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ABSTRACT 

A  distributed  parameter  maximum  principle  provides  the  basis  for  design 
of  optimal  lossless  impedance  ntatching  structures.  Given  an  arbitrary 
(wide-band)  source  at  r  =  -f  and  an  arbitrary  load  at  z  =  0,  the  optimiza¬ 
tion  procedure  synthesizes  a  uniform  cross  section,  dielectric  loaded 
waveguide  or  transmission  line  of  length  t  which  maximizes  the  real  power 
delivered  to  the  load. 

This  optimization  is  performed  subject  to  constraints  on  the  filler  material; 
p(z)  s  p^and  €(z)  is  in  [cj  .Cj*  •  •  •  T'**®  solutionis  easily  implemented, 

as  the  filler  consists  of  successive  intervals  of  readily  available  dielectric 
materials. 

An  iterative  numerical  procedure  for  solving  the  (split -boundary  function) 
necessary  conditions  is  given,  and  a  number  of  typical  solutions  for  trans¬ 
mission  line  and  waveguide  couplers  are  presented.  In  the  examples  stud¬ 
ied,  the  performance  obtained  using  only  within  0.2 

percent  of  that  obtained  using  a  continuum  of  dielectric  constants 
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SUMMARY 


This  research  is  based  on  the  upplication  of  the  results  of  the  maximum  principle  to  the  design 
of  optimal  transmission  lino  and  waveguide  structures. 

The  systems  considered  are  a  lossless  TEM  (transmission  line)  structure  and  a  lossless  TE 
(waveguide)  structure.  The  solution  sought  for  the  transmission  line  is  the  optimal  distribu¬ 
tion  along  the  line  for  the  capacity  per  unit  length,  C(z),  ze  (-f.O),  when  the  inductance  L(z)  = 
is  constant.  This  solution  may  be  implemented  by  filling  a  uniform  TEM  structure  with  a  die¬ 
lectric  material  having  a  varying  c(z).  In  the  waveguide,  we  seek  the  optimal  dielectric  filler 
material,  i.e,,  the  permittivity  e(z)  for  p(z)  = 

The  cost  functional,  to  be  minimized  in  this  process  of  optimal  design,  measures  how  effective 
a  "match"  the  electromagnetic  structure  provides  between  an  arbitrary  (wide-band)  source  at 
z  =  — f  and  an  arbitrary  load  at  z  =  0.  This  cost  is  the  total  power  reflected  due  to  impedance 
mismatch;  since  the  structures  arc  lossless,  minimizing  this  reflected  power  maximizes  the 
power  delivered  to  the  load. 

This  minimization  process  is  performed  subject  to  the  constraint  that  C(z)  e  {C^,  C^,  .  .  .  C^},  or 
€(z)  e  {e^,  This  constraint  insures  that  the  solution  is  practical  to  implement,  as  it 

only  requires  successive  intervals  of  dielectric  filler  material  having  readily  available  dielec¬ 
tric  constants. 

Necessary  conditions  for  optimality  of  the  rolution  (C(z)  or  e(z)lare  stated  in  the  form  of  a  dis¬ 
tributed  parameter  maximum  principle.  (The  systems  of  interest  are  distributed  over  both  spa¬ 
tial  aud  frequency  domains.)  Application  of  this  principle  to  the  TEM  or  TE  system  yields  a 
pair  of  coupled  partial  differential  equations  with  split  boundary-function  conditions.  An  iter¬ 
ative  numerical  procedure  for  solving  these  equations  is  given,  and  a  number  of  typical  solutions 
for  transmission  line  and  waveguide  couplers  are  presented. 

A  practical  design  example  was  selected  to  demonstrate  the  use  of  this  technique:  an  impedance 
match  is  desired  between  narrow  slots  forming  a  phased  array  receiving  antenna  and  stripline 
receivers  inserted  in  the  waveguides  which  connect  to  the  slots. 

The  small  height  of  the  waveguides  and  their  slot  openings  (chosen  to  fit  the  stripline)  results  in 
a  severe  mismatch  between  the  characteristic  impedance  of  the  wave  guided  in  the  waveguide 
and  the  free  space  wave  beyond  the  surface  of  the  phased  array;  the  V'SWR  is  greater  than  20:1 
across  the  band  atboresight  for  the  slot  pattern  under  consideration,  leaving  room  for  over  7  db 
of  improvement  by  effective  impedance  matching. 

The  method  used  in  this  investigation  is  an  iterative,  numerical  technique  based  on  the  Maxi¬ 
mum  Principle  of  Pontryagin,  and  coded  in  FORTRAN  for  execution  on  a  digital  computer. 

The  results  were  as  follows:  filling  the  waveguide  as  indicated  with  successive  intervals  of 
quartz  (e  =  3.78  e^)  and  alumina  (c  =  9  e^)  over  a  distance  of  1.15  inches,  reduced  the  VSWR 
(power-averaged  across  the  band)  from  20:1  to  5.4:1,  reducing  reflection  losses  from  7.5  to 
2.7  db. 
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When  the  available  dielectrics  range  between  ®nd  the  reflected  power  cost  showed 

negligible  improvement  with  the  addition  of  more  available  intermediate  dielectric  values:  the 
performance  using  only  *MAX^  within  0.2  percent  of  that  obtained  using  a  continuum 

of  dielectric  constants 

This  iterative  solution  required  about  10  minutes  of  run  time  on  an  IBM  360/65,  showing  very 
poor  convergence  from  a  simple-minded  initial  guess  because  of  the  Newton-Raphson.  iterative 
technique  employed.  Attempts  to  obtain  a  better  match  using  higher-c  materials  and  a  longer 
length  of  waveguide  failed  to  converge  in  a  reasonable  length  of  time. 

Changing  from  a  Newton-Raphson  iterative  technique  to  a  steepest  descent  gra  'lent  method  will 
considerably  improve  the  convergence  properties  of  the  iteration.  Once  the  procedure  is  ap¬ 
propriately  recoded,  it  is  predicted  that  a  fairly  good  impedance  match  solution  (VSWR  <1.5) 
can  be  obtained  using  reasonable  lengths  of  waveguide  (»3  inches)  and  reasonable  dielectric  ma¬ 
terials  (e.g.,  €  ~  16  €^). 
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DESIGN  OF  TRANSMISSION  LINES  AND  WAVEGUIDE 
STRUCTURES  USING  THE  MAXIMUM  PRINCIPLE 


I.  INTRODUCTION 

A.  PROBLEM  STATEMENT 

The  goal  of  this  research  is  to  apply  the  results  of  the  maximum  principle  of  Pontryagin 
to  the  design  of  optimal  lossless  transmission  line  and  waveguide  coupling  structures. 

The  systems  to  be  considered  are  a  lossless  TEM  (transverse  electro-magnetic,  i.e.,  trans¬ 
mission  line)  structure  and  a  lossless  TE  (transverse  electric,  i.e.,  waveguide)  structure.  The 
solution  sought  for  the  transmission  line  is  the  optimal  distribution  along  the  line  for  the  induct¬ 
ance  L(z)  and  capacity  C(z)  per  unit  length.  In  the  waveguide,  we  seelc  the  optimal  choice  of 
materials  filling  the  waveguide,  i.e.,  the  permeability  p(z)  and  the  permittivity  c(z)  for  zc  [— ^,  0]. 

The  results  of  the  maximum  principle  are  generally  stated  in  the  context  of  an  optimal  con¬ 
trol  problem  where  time  is  the  independent  variable  and  an  optimal  control  vector  u  is  sought 
(as  a  function  of  time)  which  minimizes  a  cost  functional  dependent  on  u(t)  and  x(t),  the  state  of 
the  controlled  system  for  te  [t^,  T]. 

We  will,  instead,  choose  a  one-dimensional  space  variable  zt  [-/,  0|  as  the  independent 
variable,  and  seek  an  optimal  parameter  vector  u  as  a  function  of  position  z  which  minimizes 
a  functional  on  the  u(z)  and  the  system  state  x(z). 

The  cost  functional  to  bo  minimized  in  each  example  will  measure  how  effective  a  "match" 
the  electromagnetic  structure  provides  between  an  arbitrary  source  at  z  =  -/  and  an  arbitrary 
load  at  z  =  0  for  a  given  signal.  Specifically,  let  the  source  be  a  transmitter  with  impedance 
Zg(w)  w.lth  spectral  power  density  S(aj)  as  the  "available"  power  density  (power  per  unit  fre¬ 
quency  which  could  be  delivered  into  a  perfectly  matched  load),  let  the  load  at  z  0  be  an  arbi¬ 
trary  impedance  Z{w)  (viz.,  an  antenna)  and  let  the  reflection  coefficient  looking  into  the  elec¬ 
tromagnetic  coupler  at  z  =  — /  be  called  p(w,  — /).  Of  course,  this  reflection  coefficient  is  a 
function  of  both  the  impedance  of  the  source  and  the  impedance  presented  by  the  coupler.  Then 
the  total  power  reflected  back  into  the  source  (not  transmitted  from  the  antenna)  is 

J  =  J  S(a;)|p(w,  -t)\^  dw  .  (1-1) 

If  we  normalize  the  source  to  1  watt  ( f  S(u)  dw  =  1),  .1  is  also  a  direct  measure  of  the  frac¬ 
tional  power  lost  (reflected  back  into  the  source)  when  a  particular  structure  is  used  to  match 
between  the  source  and  an  arbitrary  load.  This  is  the  functional  on  the  l,(z),  C(z),  [u(z),  e(z)| 
parameters  which  we  wish  to  minimize  by  finding  the  optimal  parameter  distribution  I.,*(z), 

C*(z)  [u*(z),  e*(z)l  for  ze  [-/,  0]. 


The  design  example  considered  during  this  research  involves  a  phased  array  receiving 
antenna  consisting  of  narrow  slots  distributed  in  a  uniform  pattern  over  a  plane  surface.  ICach 
slot  forms  the  mouth  of  a  waveguide  which  terminates  in  a  coupling  to  a  stripline  having  the  same 
height  and  width  as  the  waveguide;  integrated  circuit  receiving  electronics  is  embedded  in  the 
stripline  receiving  structure. 

This  antenna  operates  over  the  band  5  to  5.5  Gc,  using  waveguide  0.125  x  0.9  inch  in  cross 
section.  The  small  height  (chosen  to  fit  tlie  stripline)  results  in  a  severe  mismatch  between 
the  characteristic  impedance  of  the  wave  guided  in  the  waveguide  and  the  free  space  wave  beyond 
the  surface  of  the  phased  array;  the  VSWK  is  greater  than  20:1  across  the  band  at  borcsight  for 
the  slot  pattern  under  consideration,  leaving  room  for  over  7  db  of  improvement  by  effective 
impedance  matching. 

The  method  used  in  this  investigation  is  an  iterative,  numerical  technique  based  on  the 
Maximum  Principle  of  Pontryagin,  and  coded  in  FOHTRAN  for  execution  on  a  digital  computer. 

The  solution  obtained  lies  in  the  choice  of  dielectric  materials  which  fill  successive  inter¬ 
vals  of  the  waveguide  in  order  to  maximize  the  total  (wide-band)  power  delivered  to  the  strip¬ 
line  receiver. 

B.  BACKGROUND 

1.  Theoretical  Limitations  on  Impedance  Matching 

A  great  deal  of  effort  has  been  devoted  to  the  problem  of  wide-band  impedance  matching  in 
the  last  20  years;  we  will  attempt  to  survey  some  of  the  most  important  results  available  in 
the  literature. 

The  problem  under  consideration  is  shown  in  Fig. 1-1:  given  an  arbitrary  load  impedance 
and  a  1-ohm  resistive  source,  we  want  to  find  a  lossless  coupling  network  which  maximizes 
the  real  power  delivered  from  the  source  into  the  load.  It  can  easily  be  shown  that  this  power 
is  maximized  when  the  impedance  Z  of  Fig.  1-1  (the  impedance  seen  "looking  into"  the  lossless 
network)  is  equal  to  1  ohm^  hence  the  expression  "matching"  impedances.  Unfortunately,  it  is 
in  general  impossible  to  obtain  the  perfect  match  Z  =  1  ohm  over  a  band  of  frequencies;  we  can 
only  make  Z  approximate  1  ohm  over  a  band,  with  perhaps  Z  exactly  equal  to  1  ohm  at  a  finite 
number  of  distinct  frequencies  within  the  band. 

SOURCE  HAVING  AVAILABLE  l5-4Mli4Tl 
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Fig.  1-1.  Bosic  impedance  matching  problem. 


A  precise  measurement  of  the  quality  of  our  ''approximate  match"  over  the  frequency  band 
is  needed:  the  quality  measurement  we  choose  is  the  fractional  power  reflected: 

1  -  P, /P  = 

L'  o 

where  P^  is  the  actual  power  delivered  into  the  load  and  P^  is  the  available  power  of  the  source, 
i.e.,  the  power  which  could  be  delivered  into  an  ideal  1-ohm  load, 
p  is  the  reflection  coefficient  defined  by 

p  =  (2.  -  1)/(Z  +  1  .  (-1  <  p  <  1)  .  (1-3) 

2 


small  over  the  frequency  band  of 
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To  minimize  the  reflected  power,  we  want  to  keep  jp 
interest. 

The  first  general  theoretical  results  for  such  a  problem  were  obtained  by  H.  W.  Bode‘'  when 
he  dealt  with  the  problem  of  impedance  matching  between  a  parallel  R-C  circuit  (load)  and  a 
reSiStive  source  by  means  of  a  lossless  coupling  network.  The  fundamental  limitation  on  the 
performance  of  such  a  coupling  network  was  shown  to  be 
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where  p  i.s  the  reflection  coefficient  corresponding  to  the  impedance  Z  in  Fig.  1-1. 

If  IpI  is  kept  cot.stant  and  equal  to  |p|  over  a  frequency  bandwidth  of  W  radians/second, 
and  is  made  unity  over  the  rest  of  the  frequency  spectrum  (total  power  reflection)  then  Bode's 
integral  constraint  reduces  to 


W  In 


1 


'  max 


RC 


i.e.,  the  fractional  power  in  the  frequency  band  lost  through  mismatch  (reflection  back  into  the 
source)  is  greater  than,  or  equal  to,  exp  -(2jr/WRC).  Clearly,  this  loss  can  be  reduced  to 


zero  when  the  bandwidth  W  is  reduced  to  zero,  but  for  a  nonzero  bandwidth  it  is  impossible  to '’j] 
obtain  a  perfect  match. 

3 

R.  M.  Fano  extended  the  analysis  of  theoretical  limitations  on  broadband  matching  to  arbi 


trary  load  impedances  specified  in  terms  of  their  complex  ix>les  and  zeros, 
gral  constraints  of  the  form: 
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P(w)  In 


1 


dw  =  F(s^,  s^.  •  •  • 


whore  P(w)  is  a  ratio  of  polynomials  in  w,  p  is  again  the  reflection  coefficient  corresponding 
to  Z  of  Fig.  1-1,  and  Sj  (i  =  1,  2,  ...  n)  are  the  complex  poles  and  zeros  of  the  load  impedance. 
These  integral  equations  express  the  conditions  of  physical  realizability  of  the  matching  net¬ 
work,  and,  like  Bode's  integral  constraint,  yield  useful  information  on  the  tolerance  and  band¬ 
width  of  possible  match. 

Fano  next  made  the  point  that  when  we  are  equally  concerned  about  power  transfer  at  every 
frequency  in  a  band  W,  and  not  at  all  concerned  about  matching  outside  this  band,  the  best 
choice  for  |p(ti;)|,  subject  to  the  integral  constraint  on  lnl/|p|,  is 
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He  obtained  inte- 
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max 


This  means  that  to  achieve  the  best  power  match  over  the  band  subject  to  the  above-mentioned 
integral  constraint,  we  want  a  uniformly  good  match  within  the  band  and  total  reflection  outside 
the  band.  This  reflection  coefficient  behavior  is  shown  in  Fig.  1-2.  Note  that  with  this  strategy. 


REFLECTION  COEFFICIENT  ^(ai)  OSETS 


TOTAL  REFLECTION 


UNIFORM  COUPLING 


Fig.  1-2.  Optimal  choice  of  |p(u)|  subject  to  In  l/|p  |  du  M. 
there  is  no  point  where  a  perfect  match  is  achieved;  a  perfect  match  at  w  =  c  would  mean 


which  is  a  very  inefficient  choice  of  a  distribution  for  p,  since 


|p(w)| 


is  strictly  bounded. 

It  is  not  very  difficult  to  confirm  that  the  choice  of  uniform  |p(u))|  for  wcW  yields  a  local 
minimum  of 


subject  to  the  constraint  that  the  previous  integral  is  bounded  (see  Appendix  A). 

Unfortunately,  if  we  want  to  synthesize  a  matching  network  with  a  small  number  of  elements 
it  is  impossible  to  obtain  the  perfectly  square  |p(w)l  characteristic  as  shown  in  Fig.  1-2;  we 
can  only  obtain  an  approximation  to  it.  The  approximation  recommended  by  Fano  is  derived 
from  the  Tchebycheff  polynomial  T  (cu/o)  ). 


This  is  the  low  pass  filter  characteristic 
shown  in  Fig.  1-3;  the  reflection  coefficient  p 
oscillates  between  two  given  values  in  the  pass 
band,  and  asymptotically  approaches  unity  in 
the  attenuation  band.  For  a  bandpass  filter, 
the  argument  would  be  replaced  by 

l(w/w^  - 

where  oi^  is  the  center  frequency.  If  the  cou¬ 
pling  network  is  to  be  realized  with  k  loss¬ 
less  lumped  {LC)  elements,  the  Tchebycheff 
polynomial  is  of  order  k;  by  increasing  k  to 
infinity  we  can  approach  the  square  |p(a))| 
example, 

\  In  , — do) 

J  fpT^I 
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Fig.  1-3.  Tchebycheff  approximatian  to  optimal 
(low  pass)  lp(u)|  characteristic. 

characteristic,  and,  as  Fano  demonstrates  in  an 


can  approach  its  theoretical  bound,  i.e.,  the  integral  inequality  constraint  can  become  an 
equality. 

Assuming  that  we  have  decided  to  realize  a  (suboptimal)  Tchebycheff  match  with  a  finite 
number  of  elements,  we  need  to  find  the  element  values  which  correspond  to  such  a  match.  In 
order  to  do  this,  the  load  is  expressed  as  a  Darlington  realization,  consisting  of  a  lossless  net¬ 
work  terminated  in  a  1-ohm  resistor.  The  coupling  network  plus  load  will  then  consist  of  two 
lossless  networks  in  cascade  terminated  in  the  resistor,  as  shown  in  Fig,  1-4.  The  transfer 
matrix  of  the  second  lossless  network  is  known,  and  the  desired  Tchebycheff  transfer  matrix 
for  the  two  cascaded  lossless  networks  is  known,  so  deterniinatio.T  of  the  required  transfer 
matrix  of  the  coupling  network  becomes  a  problem  in  matrix  inversion.  Now  that  we  have  the 
transfer  matrix  of  the  coupling  network,  standard  synthesis  metnods  lead  us  to  the  element 
values. 

4 

This  approach  is  examined  in  great  detail  by  Schoeffler  who  establishes  sets  of  conditions 
under  which  the  matching  network  obtained  in  this  fashion  wil.'  be  physically  realizable,  and 
conditions  for  uniqueness  of  the  realization. 
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Fig.  1-4.  Darlingtan  realizatian  af  Tchebycheff  transfer. 
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2.  Practical  Approximations  to  Lossless  Optimal  Match 


Various  approximations  to  the  optimal  square  |p{a,-)|  characteristic  have  been  used  in  prac¬ 
tice.  For  example,  Vassiliades^  has  achieved  the  matching  network  solution  for  a  Taylor- 
Butterworth  and  for  a  Tchebycheff  characteristic  by  digital  computer  solution  of  the  algebraic 
equations  resulting  when  the  load  is  expressible  as  a  Darlington  ladder  network.  Ilis  results 
yield  the  element  for  a  lossless  matching  network  in  ladder  form. 

Other  work  where  the  lumped-coupling-element  values  are  found  for  a  Tchebycheff  transfer 

characteristic,  usually  involving  digital  computer  solution  of  the  algebraic  equations  in  con- 
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tinued  fraction  (Cauer)  form,  includes  the  work  of  Dan  Varon,  Matthaei  and  Kinariwalla. 

It  is  not  necessary  that  the  lossless  coupling  network  for  a  Tchebycheff  transfer  charac- 
teristie  be  realized  in  the  form  of  lumped  elements.  In  fact,  Dawirs  and  Tuloss  obtain  the 
desired  transfer  with  a  matching  network  composed  of  a  uniform  series  of  quarter-wave  segments 
of  transmission  line,  with  their  junctions  shunted  by  adjustable  "stubs,"  i.e.,  segments  of  shorted 
transmission  line. 

Young^°  achieves  a  Tchebycheff  characteristic  by  using  a  single  quarter-wave  segment  of 
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inhomogeneous  (different  width)  waveguide  to  couple  between  two  other  waveguides.  Cohn 


obtains  this  transfer  by  a  series  of  quarter-wave  segments  of  transmission  lines  with  different 


0  w 

Fig.  1-5.  Jacobian  elliptic  function  approximotion 
ta  optimal  |p(u)  |  characteristic. 


characteristic  impedances.  When  the  source 

R  and  load  R,  are  resistive  with  R  <  R,  , 
o  D  o  Iv 

the  characteristic  impedances  of  the  succes¬ 
sive  segments  form  a  monotonic  sequence  with 

«o<Zl<Z2-<Zn<»L 

th 

wliere  is  the  impedance  of  the  k^*^  segment 
of  transmission  line,  counting  from  the  source. 
If  R^  >  Rj  ,  the  sequence  is  monotonic  in  the 
other  direction. 

Some  work  has  been  done  using  elliptic 
functions  to  approximate  the  are  |p((ij)| 


characteristic.  The  Jacobian  elliptic  function  of  modulus  k,  Sn(w,  k)  yields  a  transfer  which  is 
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equiripple  both  in  and  outside  of  the  passband,  as  shown  in  Fig.  1-5.  Fano  discusses  this  ap- 
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proximation  with  respect  to  lumped-element  lossless  matching  networks,  and  Matthaei 


demonstrates  that  this  same  transfer  may  be  realized  by  a  nonuniform  sequence  of  quarter-wave 


transmission  line  segments. 


3.  Approximations  to  Optimal  Refiectionless  Matching 

The  previous  matching  networks  have  been  approximations  to  the  |p(u.')!  transfer  which  is 
uniform  in  the  passband  and  unity  (total  reflection)  outside  the  passband.  If  we  now  relax  the 
specification  that  the  matching  network  be  lossless,  wo  can  obtain  the  (generally  desirable) 
feature  of  no  reflection  from  the  matching  network  at  any  frequency.  This  situation  is  shown 
in  Fig.  1-6,  with  the  Darlington  realization  of  the  lossy  matching  network.  (A  lossless  network 
terminated  in  a  1-ohm  resistor.)  The  impedance  seen  looking  into  Z  is  now  1  ohm  at  all  fre¬ 
quencies,  so  the  source  is  perfectly  matched  and  no  power  is  reflected  back  into  the  source. 
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Fig.  1-6.  Lossy  reflectionless  matching  circuit. 

The  power  transfer  into  tlie  load  still  nppro.ximates  the  .square  eharaeteristie,  .so, 

by  conservation  of  jrower  the  power  outside  the  passband  and  the  power  reflected  "from  the  load 
in  the  passband  are  dissipated  in  the  "toss  resistor"  of  the  lossy  coupling  network.  The  matching 
network  may  thus  be  viewed  as  a  crossover  network  delivering  mo.st  of  the  input  power  into  the 
load  resistor  or  the  "loss"  resistor,  depr  nding  on  the  input  frequency. 

Carlin  and  Larosa^^’^^  have  dealt  with  the  necessary  synthesis  procedures  in  tletail, 
achieving  a  procedure  which  is  an  extension  of  thost'  used  to  design  the  lossless  matching  net¬ 
works  described  previously;  their  procedure  draws  heavily  on  duality  to  obtain  conjugate  im¬ 
pedances,  Their  results  include  the  fact  that  the  overall  power  transfer  into  the  load  using  a 
lossy  reflectionless  matching  network  is  at  worst  3  db  less  than  that  which  can  be  obtained  using 
a  lossless  matching  network.  In  other  words,  we  pay  a  price  of  up  to  3  db  in  power  delivered 
into  the  load  to  obtain  the  nice  feature  of  a  perfect  1-ohm  input  impedance,  Ligomenides' ^  also 
dealt  with  this  reflectionless  matching  technique,  giving  applications  in  amplifiers  and  antenna 
systems, 
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Walter  Ku  compared  the  theoretical  limitations  of  lossle.ss  matching  and  lossy  reflection- 
less  matching,  and  showed  that  the  power  transfer  of  the  lossy  network  is  exactly  3  db  less  than 
that  of  its  lossless  equivalent.  In  other  words,  to  obtain  reflectionless  matching,  half  the  power 
which  could  have  been  delivered  into  the  load  must  b<'  dissipated  in  the  "loss"  resistor.  This 
may  be  too  high  a  price  to  pay,  depending  on  the  application. 

4.  Generalization  to  an  Empirically  Measured  Load 

All  the  synthesis  techniques  we  have  examint'd  depend  on  having  the  load  impedance  ex¬ 
pressed  in  terms  of  its  complex  poles  and  zeros  at  the  outset.  However,  particularly  in  the 
ease  of  antennas,  this  information  is  not  known,  and,  in  fact,  there  is  no  reason  to  suspect  that 
the  antenna  impedance  is  a  rational  function  of  complex  frequency  s.  In  practice,  the  antenna 
imp('danee  can  be  measured  at  a  finite  number  of  frequencies 

s  --  ju-  j,  ja'2,  .  .  .  ju'^  . 

To  utilize  any  of  thi‘  i)revious  methods,  a  ratio  of  polynomials  in  s  must  be  found  which  agrees 
with  the  n  data  points,  and  this  function  used  as  a  model  for  the  antenna. 
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.Sulzer,  howevi  r,  eliminates  this  modeling  procedure  by  developing  an  algorithm  which 
yields  the  two-  or  three -lumped-element  lossless  matching  network,  based  on 
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which  minimizes 

n 

'''  1  - )  —  i  ohm  I  ^ 

il 

This  cost  functional,  magnitude-squared  impedance  deviation  from  match,  is  not  equal  to 
reflected  power,  but  is  a  first-order  approximation  to  reflected  power  when  Z  ~  1  ohm.  This 
generalization  to  an  empirically-measured  load  is  significant,  but  unfortunately  the  matching 
network  is  restricted  to  be  either  two  or  three  lumped  elements. 

5.  Generalization  to  Mixed  Distributed -Lumped  Systems 

The  matching  networks  described  previously  are  composed  either  of  lumped  elements  or  of 

quarter-wave  segments  of  distributed  structures  (transmission  line,  waveguide).  A  number  of 
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other  coupling  techniques  are  discussed  in  a  handbook  by  Matthaei,  Young  and  Jones,  including 
Tchcl)ycheff  and  Hutterworth  tr.ansfcrs  obtained  with  a  sequence  of  half-wave  transmission  line 
segments. 

Xnother  interesting  technique  described  uses  both  distributed  structures  and  lumped  ele¬ 
ments  in  tlie  coupling  network  to  obtain  the  desired  transfer:  the  structure  consists  of  a  series 
of  segments  of  uniform  transmission  line,  with  either  series-lumped  capacitors  or  parallel- 
lumped  inductances  at  each  junction.  The  capacitors  (inductances)  have  different  values  and 
the  transmission  line  segments  are  of  varying  lengths,  usually  somewhat  shorter  than  a  half¬ 
wavelength  for  the  examples  given  in  the  text.  The  appropriate  lumped -element  values  and 
lengths  of  the  segments  are  found  algebraically  from  the  coefficients  of  the  Tchebycheff 
(Hutterworth,  etc.)  polynomial.  Of  course,  a  dual  realization  could  be  used  instead,  with  paral¬ 
lel  lumped  capacitors  or  series  lumped  inductances. 

'I'his  technique  is  extremely  appealing  for  matching  into  antennas  in  the  20  to  200  Me  fre¬ 
quency  range:  the  match  can  b('  achieved  merely  by  connecting  series  or  parallel  capacitors 
to  the  leail-in  coaxial  cable  at  appropriate  points:  the  overall  structure  is  simple  and  reliable. 

6.  Graphical -Geometric -Algebraic  Techniques 

A  groat  number  of  references  appear  in  the  literature  describing  impedance  matching  tech¬ 
niques  that  use  graphical  aids  or  methods  of  analytic  geometry  (Hefs.  20  —  24,  for  example) 

These  approaches  provide  very  simple  means  of  synthesizing  an  impedance-matching  net¬ 
work  which  provides  a  perfect  match  at  a  single  frequency,  or  (when  the  load  is  purely  resistive) 
a  perfect  match  at  a  number  of  discrete  frequencies.  However,  we  have  seen  that  for  maximum 
power  transfer  wo  generally  do  not  want  to  obtain  a  perfect  match  at  any  frequency.  Further, 
we  are  concerned  about  the  match  oljtained  over  a  continuous  band  of  frequencies  while  the 
graphical  and  geometric  techniques  are  appropriate  for  analyzing  response  for  only  a  single 
frequency  at  one  time,  or  at  best,  a  number  of  discrete  frequencies  simultaneously. 

7.  Techniques  Using  a  State  Variable  Approach 

Several  groups  have  used  state  variables  in  the  design  of  transmission  lines.  For  example, 

Z  S  Z£) 

Hohrer,  Hesh  and  Hoyt  ’  investigate  the  problem  of  designing  a  los.sy  transmission  line  for 
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waveform  generation  or  pulse  shaping;  an  input  ( r(’sistivi')  signal  source  is  specified  at  z  -  — f 
and  the  wav~form  delivered  into  a  resistive  load  at  z  r  0  is  to  best  match  (in  an  integrated  quad¬ 
ratic  sense)  the  desired  output  waveform.  'I'lu'  l{(z),  l.(z),  C  (z)  and  G'z)  arc  each  constrained 
to  lie  in  an  interval;  sometimes  hang-bang  solutions  are  found,  but  usually  singular  trajectories 
occur,  whco  e  the  parameters  take  on  values  inside  tiieir  allowalile  interval. 

This  proc('  lure  was  coded  as  a  steepest  descent  iterative  procedure  in  parameter  space. 

Moth  because  of  the  singular  solutions  and  the  necessity  for  numerical  integration  of  the  state 
and  co-state  p.ai’tial  differential  equations,  convergence  was  reported  to  be  quite  slow,  typically 

requiring  5-1/2  hours  per  problem  on  the  Control  Data  1604  computer. 

27  28 

Moyer,  Wohlers  and  Kopp  ’  use  a  state  varial)le  approach  to  design  transmission  line 
structui  es  useful  for  impedance  matching  or  filtering,  selecting  the  real  and  imaginary  parts 
of  the  line  impedance  as  state  variables.  I’nfortunately,  this  choice  of  state  variables  leads  to 
a  nonlinear  set  of  partial  differential  equations  for  the  state  vector,  requiring  numerical  inte¬ 
gration.  The  parameters  Ij,  C,  and  G  (the  design  functions)  are  each  constrained  to  lie  in 
some  interval.  As  in  the  study  descriljed  above,  solutions  .sometimes  have  a  bang-bang  char¬ 
acter,  but  usually  correspond  to  singular  extremals.  This  procedure  also  utilizes  an  iterative 
gradient  or  steepest  descent  in  the  parameter  function  space. 

Considerable  difficulty  is  reported  with  convergence  to  the  singular  exti'emal;  typical  run 
times  arc  5  to  1  5  minutes  on  an  IMM  7094,  with  the  convergence  to  the  final  solutions  still  some¬ 
what  in  doubt  at  the  end  of  that  time.  Mfficient  numerical  computation  of  singular  extremals  is 
still  considered  to  be  a  prol)lcm  for  further  research. 

C.  RELATION  OF  PRESENT  WORK  TO  PREVIOUS  WORK 

Most  of  the  work  in  tlie  literature  surveyed  in  Sec.l-M  .suffers  from  several  shortcomings: 
first,  only  the  techniques  in  llcfs.  19  and  25  to  28  permit  direct  use  of  experimentally  obtained 
source  and  load  impedance  values  in  tabular  form.  All  other  techniques  require  finding  some 
model  for  the  load  tliat  would  fit  the  measured  data  points.  Unfortunately,  the  technique  of 
Ref.  19  constrains  the  matching  network  to  be  only  three  lumjaed  elements,  which  is  severely 
restrictive. 

Secondly,  only  the  technique  of  Refs.  27  and  28  actually  minimizes  the  reflected  power  of  a 
given  signal:  I'ano's  square  |p(a.')!  characteristic  is  optimal  when  the  signal  to  bo  transmitted 
has  a  uniform  power  density  for  acW  and  zero  power  for  a/W,  but  this  is  a  rather  restrictive 
class  of  signals  to  be  transmitted.  Even  in  that  case,  the  synthesis  procedures  only  approximate 
the  square  lp(a-)|  transfer  by  Tchebycheff  Taylor-Butterworth,  or  S^(u.',  K)  characteristics: 
these  choices  seem  rather  arbitrary,  and  there  is  no  reason  to  believe  that  any  of  them  actually 
m.inimizes  the  reflected  power  subject  to  the  appropriate  constraints. 

However,  even  the  technique  of  Refs.  27  and  28  has  its  drawbacks:  the  transmission  line 
parameters  L(z)  and  C(z)  are  only  constrained  to  lie  in  some  interval,  so  the  resulting  solution 
generally  has  a  continuously  varying  L(z)  and  C(z).  Such  a  structure  is  quite  difficult  to  con¬ 
struct.  A  second  difficulty  with  this  procedure  is  computational  complexity,  as  one  is  required 
to  numerically  integrate  nonlinear  partial  differential  equations. 

The  present  work  eliminates  these  drawbacks  by  placing  constraints  on  L(z)  and  C(z)  so 
that  the  solutions  are  piecewise  constant;  the  structure  can  then  be  synthesized  easily  by  placing 
intervals  of  various  commonly  available  dielectric  materials  between  the  transmission  line 


conductors.  In  addition,  a  clever  choice  of  state  variables  (in  combination  with  the  constraints 
mentioned  above)  results  in  constant-coefficient  linear  differential  equations  which  may  be  solved 
using  transition  matrices,  eliminating  all  necessity  for  numerical  integration. 

Third,  the  problem  is  generalized  to  include  TM  (waveguide  structures)  in  addition  to  the 
TKM  (transmission  line  structures)  considered  previously. 

However,  the  proi'cdure  Hi  merely  a  numerical  procedure,  so,  in  general,  the  burden  is 
on  the  user  to  decide  when  the  cost  and  the  solution  (parameter  function)  have  converged  suf¬ 
ficiently  for  his  requirements.  In  addition,  the  user  is  required  to  modify  his  starting  param¬ 
eters  in  eases  where  the  procedure  fails  to  converge, 

D.  OUTLINE  OF  PRESENT  WORK 

The  approach  taken  to  this  problem  is  the  following:  since  the  systems  being  considered 
are  distributed  through  both  a  spatial  domain  and  a  froqu('ncy  domain,  we  must  use  some  version 
of  a  distributed  parameter  maximum  principle. 
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In  See.  II,  such  a  distributed  parameter  maximum  probletn  is  developed  (following  Wang  ) 
which  provides  sufficient  conditions  for  optimality.  In  the  second  part  of  Sec.  II,  a  simple  var¬ 
iational  argument  yields  a  set  of  nccessar>  conditions  for  optimality.  These  two  sets  of  condi¬ 
tions  arc  the  basis  of  all  later  work. 

It  is  possible  to  make  a  "transmission  line  equivalent"  for  the  significant  aspects  of  the 
waveguide  behavior,  so  mathematical  descriptions  of  the  two  that  have  precisely  the  same  form 
can  be  provided.  In  faet,  with  a  proper  choice  of  state  variables  in  the  two  eases,  each  may  be 
described  by  a  linear  partial  differential  equation  and  linear  boundary  conditions  on  the  slate. 
^'urthor,  the  matrix  multiplying  the  state  in  the  partial  differential  equation  is  linear  in  the  con¬ 
trol  variable  in  each  ease.  I’or  this  reason,  Sec.  Ill  develops  the  implications  of  the  necessary 
conditions  of  See.  11  to  such  a  general  system  as  deseribed  above;  these  results  can  then  be  par¬ 
ticularized  to  the  Tl'.’M  or  Tl')  case  by  an  appropriate  sulistitution  of  parameters. 

In  See.lN',  the  lossless  TKM  structur<'  is  considered  and  a  set  of  state  variables  selected 
which  meets  the  conditions  of  linearity  describt'd  in  Sec.  111.  The  specific  line  parameters  are 
substituted  into  the  results  from  .Sec.  Ill,  and  a  set  of  necessary  conditions  for  optimality  are 
given  for  the  transmission  line. 

A  parallel  treatment  is  provided  in  See.  V,  yielding  necf'ssary  conditions  for  optimality  in 
the  waveguide  problem. 

Section  \T  provides  a  description  of  the  digital  computer  program  written  to  mechanize  an 
iterative  numerical  solution  to  th<'  necessary  conditions  given  in  Sees.  IV  and  V.  Section  VII 
contains  a  discussion  of  the  limitations  of  this  iterative  procedure,  in  particular,  suboptimality 
of  solutions  and  possible  failure  of  convergence,  and  provides  suggestions  for  minimizing  the 
effects  of  these  limitations. 

Section  \T1I  provides  the  results  obtained  by  actual  computer  runs  on  an  S!)S-940  and  an 
lHM-360  computer  for  designing  both  transmission  line  and  waveguide  coupling  structures. 
Included  are  int('r]n-etations  of  tln'  results,  conclusions,  anil  suggestions  for  improving  the  nu¬ 
merical  procedure.  Suggestions  for  future  research  are  givc-n  in  Sec.  IX. 
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II.  OPTIMAL  CONTROL  OF  DISTRIBUTED  PARAMETER  SYSTEMS 


A.  INTRODUCTION 

The  systems  considered  arc  transmission  line  (TEM)  and  waveguide  (TE)  structures.  The 
behavior  of  such  systems  is  characterized  over  an  interval  in  two  independent  variables:  space 
and  frequency.  The  space  interval  is  the  physical  length  of  the  structure  zc  [-f,  0],  and  the  fre¬ 
quency  interval  is  that  band  of  complex  frequencies  weil,  in  which  we  arc  attempting  to  propagate 
real  power  through  the  coupling  structure. 

Since  there  arc  two  independent  variables,  the  behavior  of  the  system  will  bo  characterized 
by  a  partial  differential  equation  describing  the  evolution  of  the  state  of  the  system.  Such  a  sys¬ 
tem  is  called  a  distributed  parameter  system,  and  optimization  of  the  distributed  parameters 
must  be  approached  through  a  distributed  parameter  maximum  principle. 

The  two  structures,  TEM  and  TE,  can  be  described  by  partial  differential  equations  of  thi' 
same  form;  in  effect,  we  are  able  to  find  a  "transmission  line  equivalent"  describing  the  es¬ 
sential  features  of  the  waveguide  behavior.  For  this  reason,  we  first  considt'r  in  general  the 
problem  of  optimization  of  distributed  parameter  systems,  obtaining  a  form  of  a  distributed 
parameter  maximum  principle.  In  later  sections  this  principle  is  applied  to  the  spc('ific  systems 
of  interest,  the  TEM  and  TE  structures. 
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Results  on  optimization  of  distributc'd  parameter  systems  are  available  in  the  literature. 
However,  the  systems  considered  arc  defined  ov('r  an  interval  of  time  and  space;  the  system 
is  v’iewcd  as  distributed  through  space  (scS)  and  evolving  through  time  tc  [t^,  T]  so  tlv  system 
is  defined  over  Sx  (t^,  T].  This  same  notation  will  be  used  in  the  following  developments  for 
notational  convenience;  however,  we  will  be  able  to  rcinti'rpret  the  final  results  by  substituting 
the  independent  variables  of  interest,  wcR  and  ze  [-i,  0)  for  s  and  t,  respectively,  obtaining  a 
system  defined  over  fix  (-i,  0].  In  effect,  u  is  th<.‘  "spatial  distribution"  variable  and  z  is  the 
independent  "time"  variable. 

B.  PROBLEM  STATEMENT 

We  consider  the  following  system: 

(1)  A  distributed  parameter  system  whose  state  at  time  t  is  given  by  x(t,s), 
a  vector-valued  function  defined  over  te  (t  ,  T]  and  scS,  a  compact  subset 
of  Euclidean  N-spacc. 

(2)  The  state  of  the  system^  evolves  through  time  according  to  the  partial  dif¬ 
ferential  equation: 


tOne  example  illustrating  this  type  of  system  is  the  heat-equation  relation  for  a  uniform  piece  of  material  with 
internal  (volume)  heat  sources: 


9x(t,s) 

at 


K  a^x  t,s  .  1  , 

pC  g2  pC 


where  x(t,s)  is  the  temperature  in  ®K,  u(t,s)  is  the  rate  of  internol  heat  generation  in  cal/cmysec,  p  Is  material 
density  in  gm/cm^,  C  is  specific  heat  in  cal/cm-®K,  and  K  is  the  thermal  conductivity  in  cal/sec-cm-'K. 
Clearly,  f  is  here  a  linear  scalar  operator  involving  a  second  space  partial  derivative  on  the  state  x  plus  a  con¬ 
stant  times  the  control  u. 
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(2-1) 


9x(t,  s) 

at 


f  (x(t,  s),  ll(t,  s).  t,  s] 


when.'  u  is  thf  control  strategy,  a  vector-valued  function  also  dc'fined 
nv(,‘r  [tg,  'I’l  xS  and  constrained  to  be  in  the  set  of  admissible  controls 
*J(lt().  (’]  xS),  and  f  is  a  given  operator  involving  algebraic  operations 
on  the  control  u  and  st)atial  partial  derivatives  on  the  state  x. 

(3)  Tliere  is  a  scalar  cost  functional  or  performance  index  on  the  control 
strategy  u([t^xT],  S)  which  may  be  represented  as 

.1  {x^(s).  t^.  u([t^.  T]  xS)}  =  £  K  (T,  s,  x(T,  s)]  dS 

^  L  |t,  s,  x(t,  s),  u(t,  s)]  dS  dt  .  (2-2) 


Tlic  first  term  represents  a  spatially  integrated  terminal  time  cost 
(a  functional  on  the  terminal  state)  and  the  second  term  represimts 
the  time  integral  of  a  cost  rate  which  is  itself  a  spatial  integral. 

Given  the  system  specified  by  (1),  (2)  and  (3),  an  initial  state  x^^(s)  =  x(t^,  s)  and  a  terminal 
timi'  T,  we  want  to  find  the  optimal  control  u(t,  s)  for  tc  [t^,  T]  and  scS  wliich  results  in  the  min¬ 
imum  value  of  tile  cost  functional: 

J’Mx^(s),tl=  min  .1  {x  (s),t  .u((t  ,TlxS))  .  (2-3) 

°  °  ueU([t^,Tl.xS)  o  o  o 


C.  HAMILTON-JACOBI  EQUATION  IN  FUNCTION  SPACEt 

Bellman's  Principle  of  Optimality  states  thtit  any  portion  of  an  optimal  trajectory  in  state 
space  (th<'  trajectory  described  by  the  given  system  when  the  optimal  control  is  applied)  is  itself 
an  optimal  trajectory  between  tlie  initial  and  terminal  states  of  that  portion.  If  this  were  not 
true,  a  Ic'ss  costly  control  would  e.xist  for  that  portion  of  the  trajectory'  and  thus  for  the  entire 
trajectory;  tills  violates  the  hypothesis  that  the  original  trajectory  is  optimal. 

Applying  this  principle  of  optimalitv  to  the  problem  stated  above,  with  the  portion  of  tra¬ 
jectory  under  consideration  being  x(t,  s)  for  scS,  te  (t,t  +  e],  wo  have: 

.1’^“  [x(t,s),  t)  =  min  j.J*  (x(t  +  €,s),  t  +  c] 

ueU  [t,  T  +  e)  xS  | 


pT  +  e  r 

+  ^  ^  L  (t,  s,  x(t,  s  ),  u(t,  s))  dS  dt 

where,  of  course,  x(t  +  e,  s)  depends  on  u(t,  s)  for  tc  (t,  t  +  e],  scS.  Thus, 

T+e  n 

^  L  [t,  s,  x(t,  s),  u(t,  s)]  dS  dt 

‘S 


[x(t,  s),  t)  -  J*  [x(t  +  €,  s),  T  +  ej  < 


(2-4) 


(2-5) 
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t  Adapted  from  Wang. 


Wo  now  make  use  of  the  following  assumptions; 

(1)  L  (t,  s,  x(t,  s),  u(t,  s)]  clS  is  bountictl  from  below. 

(2)  f  [t,  s,  ,x(t,  s),  u(t,  s)]  is  a  sufficiently  smooth  operator  so  that  for  a  given 
initial  state  and  control,  the  system  partial  differential  Eq.  (2-1)  has  a 
unique  solution. 

(3)  J*  [x(t  s),t]  is  smooth  in  t,  i.e.,  obeys  the  Lipschitz  condition 
I  9J  ^9t1  <  ^  . 

(•1)  J*  (x(t,s),t]  is  a  smooth  functional  on  x(t,s). 

Under  these  assumptions,  the  equation  above  may  be  expanded 


-  ^  [x(t,s).t]4  <[  —  .  f  [X(T,S),U(T,S),T,S]>J 


rj(x  +  hj)- J(x)l 

lim  — -T  — -  (row  vector)  (2-7) 

Jg  |h.(s)|dS-o^  Jg  hj(s)  d.S 

wliere  h.(s)  is  any  continuous  vector  function  (of  the  same  dimension  as  x)  whose  only  nonzero 
component  is  the  i^^  component,  and  which  lias  compact  support  in  a  region  AS  in  S,  (i.e., 
h(s)  =  0  for  s^AS,  h(s)  ^  0  for  seAS).  The  requirement  tliat  J  [x(s,  t  ),  t)  be  smooth  in  x(s,  t  ) 
is  necessary  for  the  existence  of  the  limit  corresponding  to  the  above  expression. 

The  llamilton-Jacobi  equation  we  have  derived  may  be  written  more  compactly  if  we  define 
the  costatc  column  vector: 

p(t,  s)=  [  ~ax'(I ’s')^*'l  (column  vector)  (2-8) 

and  the  Hamiltonian 


fi.I  (x(t,  s ),  t] 
'5.x(t,s) 


Il(t,  p,  x,u)=  <p(t,  s),  f  [t,  s,  x(t,  s),u(t,  s))>g  +  ^ 


L  (t,  s,  x(t,  s),  u(t,  s)]  dS 


wherein  the  inner  product  brackets  again  denote  the  integral  over  S  of  the  vector  inner  product 
of  p  and  f 


Then  tlie  llainilton-.laeobi  equation  is; 


3.1*  s),tl 
5t 


ntin  ll(t,  p,  X,  u) 
utl  (txS) 


(2-10) 


with  houndaev  condition 


J*  Ix(T.s).  T1 


K  [x(  l',  s),  I'l  dS 


(2-11) 


It  is  simple  to  verify  tliat 

3x*  (t,  s)  .611*  (t,  p*  ,  X’  )  ' 

at  ■  “  ^  6p'  ■  ' 

and 

3p*  (t,  s)  _  _  I  611*  (t,  p-' ,  X*  ) , ' 

3 1  dx* 

wliorc  II*  is  the  u-minimal  II,  i.e., 

II*  (t,  p*  ,  X*  )  =  min  ll(t,  p*  ,  x*  ,  u) 
ucL  (txS) 

and  X*  and  p*  are  the  corresponding  optimal  state  and  costate  functions.  Since 


(2-12) 


(2-13) 


(2-14) 


(2-15) 


then 

(2-16) 

(assuming  that  the  indicated  Frechet  derivative  exists)  and  we  are  also  given  the  initial  state 

x(t^,  s)  =  x^(s)  .  (2-17) 

Thus  we  have  a  pair  of  partial  differential  equations  with  boundary  conditions,  and  the  op¬ 
timal  control  problem  reduces  to  that  of  solving  a  two-point  boundary  value  problem  in  function, 
space. 

These  conditions  are  sufficient  for  optimality,  i.e.,  if  the  "optimal  return  function"  J*  (x(s),tl 
is  found  which  provides  the  functions  u*  ,  x*  and  p*  satisfying  all  the  conditions  above,  any  other 
choice  of  u^u*  will  result  in  a  cost  at  least  as  great. 

D,  SUMMARY 

Sufficient  conditions  for  optimality  have  been  derived  based  on  the  following  assumptions; 


P*  (T,  s) 


6 

6x* 


K  |T.s.x*(T.  s)] 


dS 


(1)  L  [t,  s,  x(t,  s),  u(t,  s)]  is  bounded  from  below.  Any  positive  definite  cost 
rate  will  trivially  meet  this  requirement. 

(2)  f  [t,  s,  x(t,  x),  u(t,  s)|  is  a  sufficiently  smooth  operator  such  that,  for  a 
given  initial  state  and  control,  the  system  partial  differential  Eq.  (2-1) 
has  a  unique  solution.  This  requirement  i.s  satisfied  in  most  systems 

of  physical  interest,  for  example,  the  heat  equation  and  the  wave  equation. 
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(3)  J*!"  [x(t,s),t]  is  smooth  in  t,  i.e.,  obeys  the  Lipschitz  condition  |8J/9t|  <  . 

This  condition  is  equivalent  to  requiring  tliat  the  Hamiltonian  remain  finite 
for  all  te  [t^,  T]. 

(4)  J*  [x{t,  s),  t]  is  a  smooth  functional  on  x(t,  s),  i.e.,  the  Frechet  derivative 

P(T.s)  =  [  axlri  exists  for  all  te  [t^,T]  . 

In  particular,  for  t  =  T, 


P(T,  s) 


_ 5 _ 

6x{T,  s) 


r  K  {x(T,  s),  T}  dS 


exists. 


Uuder  tliesc  four  assumptions,  sufficient  conditions  for  a  globally  optimal  u’>  (t,  s)  for  t  [t^,  T] 
and  seS  are  the  following.  If  we  define 


Il{t,  p,  X.  u)  =  <p{t.  s),  f  [t,  s,  x(t,  s).  u(t,  +  \  F  (t,  s,  x(t,  s),  u(t,  s)]  dS 


(2-18) 


[x{t.  .s).  t] 


min  Ii(t,  p,  X,  u) 
uili{txS) 


a  nd 


.1*  [x(T,  s),T]  =  K  [x(T,s).TJ  dS 
or  equivali'ntly,  under  all  the  assumptions  made, 

ax*{t,^)  .  aiPMt.  p,  X)  ' 

5t  ~  ^  6p  ‘ 

§12l(l._s)  -  I  61 1’’’’  (t,  p,  x;,' 

"^t  -  I  6x  * 


(2-19) 


(2-20) 


(2-2i) 

(2-22) 


then 


ll(.x’!‘ .  p’i"  ,  u<‘  ,  t)  <  Il(x*  ,  p'’  ,  u,  t)  V  ucL’(t.xS)  ,  V  te  (t^^,  T) 


(2-23) 


with  the  equality  holding  only  for  u  =  u’'''  . 

x(tg,  s)  =  .x^(s)  (2-24) 

and 

(2-25) 

E.  NECESSARY  CONDITIONS  VIA  VARIATIONAL  METHODS 

Given  the  system  described  in  .Sec.  I,  we  can  quite  easily  state  a  necessary  condition  which 
an  optimal  control  strategy  must  satisfy;  the  cost  must  not  decrease  when  any  small  allowable 
perturbation  of  the  cost  sti'ategy  is  performed. 

The  candidate  for  optimality  is  u  (t,  s)  cU([t^,  T]  .xS);  the  iierturbed  control  is  u  (t,  s)  =  u  (t,  s) 
for  tc  [t,  t  +  e]  and  u  (t,  s)  =  w(t,  s)  for  te  [t,  t  +  c],  (wcf).  Tiie  initial  conditions  on  the  state 
[x(t^,  s)  =  x^(s)]  remain  unchanged;  since  thi;  target  set  for  t  =  T  is  the  entire  function  space, 


p(T,  s) 


K  [x(t,s),  T]  d.S 
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this  perturbed  control  satisfies  the  boundary  conditions  and  is  thus  a  candidate  for  optimality. 
The  cost 

J  {x^(s),  t^,  u([t^,  T]  xS)}  =  ^  K  [x(T,  s)l  dS  +  J  H(t,  p,  x,u)dt 

‘'o 

-  ^  p'(t,s)  (t,s)  dS  dt  (2-26) 


€ 


We  now  make  use  of  the  following  assumptions: 

»)  /s  Mt  ,  s,  X,  u)  dS  is  bounded  from  below. 

(2)  f(t,  s,  X,  u)  is  sufficiently  smooth  that  for  a  given  initial  state  and  control, 
there  Is  a  unique  solution  to  the  system  Eq.  (2-1). 

(3)  The  Hamiltonian  H(t,  x,  p,  u)  is  finite  for  finite  norm  u. 

(4)  K  (x(T,  s),  T]  dS  and  L  [t,  s,  x(t,  s),  u(t,  s)]  dS  are  continuously  dif¬ 
ferentiable  in  X,  i.e.,  the  Frechet  derivative  6/6x  of  each  exists. 

(5)  p(t,  s)  is  continuous  in  t. 

Then  by  forming  the  difference  in  costs 

x(T,  s)) 

/s 

+  ^  <[|^  (X,  p.u.t)  +  ,  6x(t,  s)  dt>s 


[Xo(s),  uj  -  J 


[X^(s) 


to.ul  = 


6x  ^ 


K{x(T,6).T}dS-p'(T) 


H(x,  p,  u,  t)  -  H(x,  p,  u,  t)  dt 


(2-27) 


To  insure  that  the  cost  difference  above  is  >0,  we  select  the  continuous  (Lagrange  multiplier) 
function  p(t,  s)  to  have  the  following  properties: 

9p'(t.  s)  _  _  (5H(x,  p,u,  t)  (2-28) 

6x 

p'(T.s)  =  r  K  [x(T.  s),T)  dS  .  (2-29) 

6x 


The  first  two  terms  are  thus  zero  for  all  permissible  u(t,  s);  the  remaining  term  must  then 
be  >0;  letting  c  -•  0  we  have: 

H(x‘' ,  p*  ,  u*  ,  t)  <  Hlx*  ,  p*  ,  u".  t)  for  all  uEU(txS)  (2-30) 


where  the  *'s  now  denote  optimal  quantities. 
The  necessary  conditions  arc  thus; 


Ox*  _  ,  (SH*  (x,  p.  u,  t),' 

~5r  - 1 — ^ — J 


(2-31) 
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(2-32) 


xM^s)  x^(s) 

9p’'  ,  6II''’ (x,  p,  u,  t),' 

ar  6^'  * 


(2-33) 


p"'  (T,  s) 


6 

6x'- 


\  K  (x’MT.  s).  T)  dS 

‘S 


(2-34) 


H(x''' ,  p”^  ,  u'''  ,  t)  <  II(x<‘ ,  p*  ,  u,  t) 

for  all  ueU(txS). 


(2-35) 


F.  DISCUSSION 

Tli(‘  difforonco  bctwocn  thosi'  nccr.ssary  conditions  and  the  sufficient  conditions  derived 
previously  is  that  in  the  suffieient  conditions  u*  uniquely  minimizes  tlic  Hamiltonian,  vvliile  in 
the  necessary  conditions,  the  minimizing  u  may  not  bo  unique.  However,  the  sufficient  con- 
dit.ions  were  derived  under  one  additional  assumption:  that  .1*  [.\(s,t),t]  is  smootli  (P'rochot 
diffen'ntiablc  in  x(s,  t)  for  all  tc  [t^,  T] } . 

If  a  system  of  interest  satisfies  all  the  assumptions  of  th(,‘  sufficient  conditions,  then  the 
solution  of  Hqs.  (2-21)  to  (2-25)  yields  tlu'  globally  optimal  control  function;  this  candidate  must 
then  be  compared  with  all  other  eandidaies  for  optimality  to  determine  which  one  is  globally 
optimal. 

There  are  two  ways  in  which  a  solution  could  moot  the  requirements  of  the  necessary  con¬ 
ditions,  but  not  those  of  the  sufficient  conditions.  First,  the  optimal  return  function  J’'‘  [x(t,  s),tl 
may  not  be  a  smooth  functional  in  x(t,  s).  In  the  finite-state  optimal  control  problems,  this  oc¬ 
curs  when  the  state-space  trajectory  strikes  the  target  set  tangentially,  i.o.,  the  trajectory 
barely  clips  the  target;  an  infinitesimal  change  in  the  trajectory  near  the  target  would  necessitate 
drastic  changes  in  the  system  behavior,  as,  for  example,  having  to  overshoot  the  target,  then 
turn  around,  and  return.  Such  a  condition  is  referred  to  as  "abnormality,"  and  means  that  the 
target  is  exactly  on  the  border  of  the  set  of  reachable  states,  and  that  an  infinitesimal  perturba¬ 
tion  in  the  initial  state  cannot  be  corrected  by  an  infinitesimal  perturbation  of  the  control  strategy. 

An  essential  feature  of  the  abnormality  condition  is  the  restricted  target  set.  In  the  dis¬ 
tributed  parameter  problem  at  hand,  we  selected  the  entire  function-space  as  the  target  set. 

It  thus  seems  plausible  that  abnormality  conditions  and  thus  "corners"  in  .1’’'  [t,  s,  x(t,  s)]  with 
respect  to  x(t,  s)  will  not  appear:  every  state  lies  within  the  target  set  so  there  is  no  danger  of 
the  trajectory  merely  "clipping"  the  target  tangentially. 

If  this  is  indeed  the  case,  we  nec'd  only  worry  about  cases  that  satisfy  the  necessary  con¬ 
ditions,  but  not  the  sufficient  conditions  because  of  the  second  reason:  the  function  u’"  (t,  s)  which 
minimizes  ll(x'‘  ,  p’i' ,  u,  t)  is  not  unique.  If  the  function  u*  (t,  s)  is  only  nonunique  over  a  sot  of 
te  [1q»  T]  with  measure  zero,  we  need  not  be  concerned,  because  any  of  the  minimizing  functions 
u  will  result  in  the  same  cost. 

If  H(x'' ,  p*  ,  u,  t)  has  some  small  number  of  minimizing  solutions  over  a  nonzero  time  inter¬ 
val,  then  a  direct  comparison  of  the  eo.sts  resulting  from  these  several  strategies  is  possible, 
and  will  yield  a  globally  optimal  strategy.  However,  if  there  is  an  infinite  number  of  minimizing 
functions  u  over  some  finite  time  interval,  we  are  faced  with  the  "singular"  condition,  and  a 
different  approach  is  needed  to  find  the  optimal  control. 
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Ono  mi'tliotl  is  to  computo  the  time  derivative  of  the  Hamiltonian 


dt 


(x*  ,  p’''  ,  t) 


an* 

at 


(x* 


3x*v  .511*  9p*.  aU 

at  ^6p*  '  ■  at 


Successive  time  derivatives  of  this  relation  will  yi<'ld  relations  between  u*  ,  x*  ,  p*  ,  and 
their  derivatives.  In  any  particular  case,  it  may  be  possible  [by  substituting  into  the  conditions 
of  i;qs.(2-31)  through  (2-35)]  to  find  a  unique  u*  which  satisfies  these  new  conditions  and  also 
satisfies  the  boundai-y  conditions. 

If  a  solution  to  the  original  problem  docs  exist,  this  is  then  the  solution. 
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UI.  LINEAR  DISTRIBUTED  PARAMETER  SYSTEMS 
WITH  LINEAR  TERMINAL  COST 


A.  INTRODUCTION 

In  the  previous  section  wc  derived  optimality  conditions  for  systems  distributed  over  space 
(scS)  and  evolving  through  time,  tc  [t^,  T),  so  the  domain  of  system  definition  is  Sx  [t^,  T].  Now 
we  wish  to  analyze  electromagnetic  structures  which  propagate  real  power  in  the  +z  direction 
(by  standard  convention),  and  thus  have  physical  extent  zr  (— <,  0).  Power  flow  will  be  considered 
for  complex  frequencies  w  in  some  frequency  interval  fi,  so  we  make  a  substitution  of  variables 
in  our  previous  results,  w  replaces  the  spatial  variable  s  and  z  replaces  the  time  variable  t. 

There  is  one  additional  point  to  be  considered,  however;  since  real  power  flow  is  in  the  +z 
direction,  we  connect  the  source  to  the  matching  structure  at  z  =  — f  and  to  the  load  at  z  =  0.  The 
cost  functional  to  be  considered  will  bo  the  power  reflected  into  the  source  at  z  =  — f.  (Since  the 
matching  structures  under  consideration  will  be  lossless,  minimizing  this  reflected  power  will 
maximize  the  power  delivered  to  the  load.)  In  the  original  derivation  the  cost  involves  a  termi¬ 
nal  cost  term  which  is  a  functional  on  the  state  at  t  T.  In  the  present  system,  the  cost  is  a 
functional  on  the  "terminal"  state  at  z  =  — f.  Wo  must  then  consider  z  =  0  as  the  initial  "time" 
and  z  =  — /  as  the  terminal  "time."  The  system  is  thus  propagating  backward  through  "time." 

This  makes  an  essential  change  in  the  necessary  conditions.  To  minimize  the  total  cost, 
the  Hamiltonian  must  be  maximized  rather  than  minimized,  as  in  Eq.  (2-35);  this  may  be  seen 
by  writing  the  cost  as  a  "time"  integral,  then  changing  the  sign  of  the  integrand  to  interchange 
the  limits  of  integration. 

As  we  noted  previously,  the  state  equations  describing  the  TEM  structure  and  those  describ¬ 
ing  the  TE  structure  have  the  same  form.  In  fact,  both  can  be  expressed  as  linear  partial  dif¬ 
ferential  equations.  In  both  cases,  the  cost  due  to  reflected  power  loss  can  be  expressed  as  a 
linear  functional  on  the  "terminal"  state  at  z  =  -f. 

Since  the  two  problems  arc  so  similar,  we  will  investigate  our  reinterpreted  necessary  con¬ 
ditions  in  the  general  case  of  a  system  with  a  linear  first-order  partial  differential  equation  on 
the  state  and  a  linear  terminal  cost.  The  results  obtained  will  then  be  applied  to  the  specific  ex¬ 
amples  of  the  TE.M  structure  and  the  TE  structure  in  later  sections. 

B.  PROBLEM  STATEMENT 

We  consider  the  following  problem: 

(1)  A  distributed  parameter  system  whose  state  at  location  z  is  given  by 
x(z,w),  a  vector-valued  function  defined  over  ze  [— f ,  0]  and  wcH  (an  in¬ 
terval  of  comple.x  frequencies). 

(2)  The  state  of  the  system  evolves  through  z  according  to  the  partial  dif¬ 
ferential  equation 

=  A  [w,  u(z)]  ■  x(z,w)  (3-1) 

where  u  is  the  control  strategy,  a  vector-valued  function  defined  over 
[— f,  0]  and  constrained  to  lie  in  the  set 

U  =  {ui,U2,...u^}  . 

I'urther,  A  [a!,u(z)]  is  a  matrix  linear  in  u  for  all  ze  [—1.  0]  and  wefl: 
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{ 


A  [ciJ,  u(z)|  =  A^(a))  +  2j  Aj(w)  u.(z) 
i=l 

All  Aj(ci))  are  continuous  in  oj.  The  initial  condition  on  the  state  is  x(0,  w)  = 
XqIw)  with  Xq{<jli)  continuous  in  o). 

(3)  There  is  a  scalar  cost  functional  or  performance  index  on  the  control  strat¬ 
egy  which  may  be  represented  as; 


J  {x  (cij),  u  [-<,  0|}  =  \  K'(aj)  x(— 1,0!)  do) 

°  •-',.,,0 


This  is  a  frequency  integrated  "terminal  location"  cost,  or  a  linear  func¬ 
tional  on  the  terminal  state,  with  K(a))  an  n-vector  continuous  in  oi. 

Given  the  system  specified  by  (1),  (2),  and  (3),  an  initial  state  Xo(u))  = 
x(0,ci))  and  a  terminal  location— f,  wc  want  to  find  the  (optimal)  control 
u(z)  for  zc  [— {,  0|  which  results  in  the  minimum  value  of  the  cost  functional: 

J  [x  (a))|  min  J  {x  (oj),  u(1-{,  0))} 

°  ueU([-<,01)  " 

C.  NECESSARY  CONDITIONS  FOR  OPTIMALITY  -  THE  HAMILTONIAN 


h’or  the  problem  outlined  above,  wc  define  the  Hamiltonian 
II(x,  p,  u,  z)  =  y  p'(z,w)  x(z,w)da) 


The  necessary  conditions  state  that  a  costate  vector  must  exist  (where  the  indicates 

an  optimal  quantity)  such  tliat: 

(1)  -  [— =  A  [a),u-!Mz)|  •  x-M^.u))  (3-6) 

with  x*(0,w)  =  Xo(w).  This  condition  merely  restates  the  partial  differ¬ 
ential  equation  of  the  system  and  its  initial  state: 

3p<MZ'_c£)  ,,  _j^.*.(x_^|.u,_z)j'  [u..u«C/.)|  •  P-'M^.O))  (3-7) 


with  boundary  condition 


p(-/,  w)  = 


. " 


'(w)  du' I'  K(a;) 


(3)  II(x’''‘,  p’' ,  u*,  z)  >  Il(x’^,  p  ,  u,  z)  forallufl  ,  i.c.,  u’> 


maximizes 


\  p'^tz.^j)'  A[(a),u(z)|'  x(z,a')da;-  \  p'’>(z,a)) 

M 

■  A  (w)  •  x’Mz,a;)da)+  ^  u.(z)  \  p'*  (z,  w)  A.(a;)  x”^' (z,  w)  do,’  .  (3-10) 


The  first  term  is  independent  of  u.  Writing  the  Hamiltonian  as 
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(3-11) 


M 

II(x>:',p’!',  u,  z)  -  +  Yj 

i  1 

this  maximization  requires; 


u. 

if 

m.(z) 

>  0 

‘max 

L 

U. 

‘min 

if 

m.(z) 

<  0 

1 

determined  only  to  lie  in  U 

if 

m.(z) 

=  0 

where  the  MAX  and  MIN  subscripts  denote  the  maximum  and  mini¬ 
mum  available  control  values  in  U.  The  third  possibility  of  Eq.  (3-12), 
mi(z)  =  0,  must  be  considered  in  more  detail.  !■  irst,  suppose  that 
mj(z)  =  0  only  over  a  set  of  z  having  measure  zero:  tlien  the  selection 
of  any  ucU  at  those  locations  will  result  in  the  same  system  trajectory 
and  total  cost.  If,  however,  nijlz)  0  over  some  interval  of  z,  the 
choice  of  u  becomes  significant.  Different  selections  of  u  will  gener¬ 
ally  result  in  different  system  trajectories  and  different  costs  because 
Eq.  (3-12)  merely  gives  a  necessary  condition  for  optimality,  not  a  suf¬ 
ficient  condition. 

This  condition,  mj(z)  =  0  over  a  nonzero  interval  of  z,  is  known  as  the 
.singular  condition.  The  appearance  of  this  singular  condidon  greatly 
complicates  any  attempt  at  numerical  solution  of  the  necessary  condi¬ 
tions,  as  we  shall  discuss  later  on. 

To  avoid  the  complexities  of  this  singular  condition  (whicJi  is  usually 
present  in  the  problems  to  bo  considered  below)  we  employ  a  mathemat¬ 
ical  artifice:  a  slight  modification  of  our  original  cost  functional,  fur¬ 
ther,  we  consider  u(z)  to  be  a  scalar  for  simplicity  in  the  following  argu 
ments.  Generalization  to  a  vector  is  straightforward,  but  unnecessary 
for  the  problems  at  hand. 


D.  MODIFIED  COST  FUNCTIONAL 
1.  Necessary  Conditions 

It  will  bo  shown  that  it  is  desirable  to  make  a  small  modification  to  the  cost  functional  of 
Eq.  (3-3).  In  particular  wo  add  a  term  such  that 

J(x  ,  0)  -  C  K'(w)  x(-f,a;)  doi  +  ^  ^  |  u  -  u  |  ^  dz  (3-13) 

The  u  term  is  some  value  of  control  intermediate  betw'oon  the  extremes  of  ucl'  =  {  u, ,  u,,  .  .  .  u  ) . 
a  '  1  ’  2’  m' 

Since  wo  are  in  reality  interested  in  the  cost  f  „  K'x  dca,  but  not  about  (u  —  u  y,  A  is  chosen 
sufficiently  small  that  the  contribution  of  the  |^|u  —  u^[  term  is  negligible  (for  engineering  pur¬ 
poses)  with  respect  to  the  /  K'xdw  term. 

Even  though  this  term  is  negligible  for  engineering  purposes,  it  is  quite  significant  mathe¬ 
matically.  The  term  is  a  mathematical  artifice  that  greatly  simplifies  consideration  of  the  so¬ 
lutions  of  the  necessary  conditions  that  correspond  to  singular  extremals;  we  would  like  to  solve 
the  problem  in  the  limit  as  A  —  0,  but  it  is  neccssai  y  to  keep  A  sufficiently  large  in  tlie  numer¬ 
ical  procedures  that  it  does  not  disappear  in  round-off  and  truncation  errors. 

If  we  rewrite  the  Hamiltonian  of  Eq.  (3-4)  to  include  this  term,  we  have 

H(x,p,u,z)=\  p'(z,cj)-  A  [w.u(z)l  ■  x(z,  w)  dw  —  A  |u  -  u  |  ^/2  (3-14) 
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and  Eq.  (3-11)  becomes 


ll(x,i),u,z)  m^  +  M  •  u  — ^|u— u,^]V2  •  (3-13) 

The  i)nrtinl  differential  equations  for  the  state  and  costatc  Eqs.  (3-6)  and  (3-7),  and  also  the 
boundary  conditions  on  the  state  and  costate  are  not  changed  by  tlie  addition  of  tlie  A  term.  It 
only  influences  the  total  cost  (by  a  very  small  amount)  and  the  Hamiltonian. 

'I’he  change  in  the  Hamiltonian  is  quite  helpful,  however:  we  note  that  the  Hamiltonian  is 
now  a  quadratic  function  of  the  control  instead  of  a  linear  function.  The  singular  condition  m^  =  0 
no  longer  prevents  the  Hamiltonian  from  having  a  well-defined  ma.ximum  with  respect  to  u. 

In  fact,  since  u  is  a  scalar,  it  is  evident  that 

u’l'  =  ur  { u, ,  u,, .  .  .  u  }  closest  to  (u  +  iM/A)  (3-16) 

'  1'  2'  m'  '  a  ' 

whore  M  is  the  factor  multiplying  u  in  the  Hamiltonian.  This  choice  of  u*  maximizes  the 
Hamiltonian  over  all  available  ucl  ,  and  this  maximization  usually  provides  a  unique  value  for  u’^'. 

However,  for  some  values  of  M  (this  sot  has  measure  zero)  two  admissible  values  of  u  are 
equally  close  to  u,^  —  M  A.  This  condition  will  now  he  called  the  singular  condition  since 
Eq.  (3-16)  docs  not  yield  a  unique  u''',  but  only  specifics  that  the  optimal  value  is  one  of  the  two. 

Note  that  if  we  had  not  used  the  artifice  of  the  A  |  u  —  u^  |  ,  Z  term,  whcnev'cr  M  =  0,  we  would 
be  faced  with  a_H  possible  controls  u^.u^.  .  .  .  u^^  as  equally  good  candidates  for  optimality. 

2.  Singular  Condition 

When 


u^  +  M/A 


“k+1 


we  are  faced  with  ambiguity  in  Eq.  (3-1 6);  both  Uj^  and  maximize  the  Hamiltonian.  If  M 
takes  on  this  value  only  over  an  interval  of  z  having  measure  zero,  then  either  control  is  equally 
good.  However,  if  M  remains  at  this  critical  value  over  some  nonzero  interval  of  z,  its  first 
and  second  derivatives  must  be  zero  in  that  interval; 


^  0  (3-17) 

dz 

2 

— ~  =  0  (3-18) 

dz 

Hecause  A  [w,u(z)]  is  linear  in  u,  Eq.  (3-18)  gives  an  e.xplicit  solution  for  u,  i.c.,  Eq.  (3-18) 
holds  if  and  only  if  u  takes  on  a  particular  value,  say  u^. 

.Several  possibilities  arise:  first,  suppose  Eq.  (3-17)  holds  and  u^  equals  cither  u^^  or 
then  a  singular  extremal  exists,  and  the  value  of  u  along  this  singular  extremal  is  Uj^  or 
as  specified  by  Eq.  (3-18).  This  resolves  the  ambiguity  of  Eq.  (3-16). 

If  either  of  the  two  above  suppositions  fails,  then  this  is  not  a  singular  extremal,  .M  main¬ 
tains  its  critical  value  only  over  a  distance  of  measure  zero,  and  it  makes  no  difference  whether 
U|^  or  is  selected  at  this  point. 

Thus,  because  of  the  added  A  I  u  —  u^  |  ^  term,  the  singular  condition  causes  no  now  difficul¬ 
ties  in  specifying  u’'‘(z)  in  the  numerical  solution  of  the  necessary  conditions  Eqs.  (3-6)  to  .(3-9). 


llowc'vcr,  as  we  will  see,  this  singular  condition  can  cause  difficulty  in  obtaining  convergence 
when  an  iterative  method  is  used  to  obtain  simultaneous  solution  of  Eqs.  (3-6)  to  (3-9). 


E.  PROPERTIES  OF  SOLUTIONS  -  THE  HAMILTONIAN 

Due  to  the  imposed  constraint  on  the  control  uc  {u^,  u^, .  .  .  any  control  function  u(z) 

which  satisfies  the  necessary  conditions  Eqs.  (3-6)  to  (3-9)  must  be  a  piecewise  constant  function 
of  z  over  [-f,  01;  it  will  switch  between  available  values  in  (  u. ,  u,, ...  u  }  over  successive  in- 
tervals  of  z.  Over  each  of  these  intervals,  the  matrLx  A  (co,u(z)]  will  be  constant  in  z. 

This  moans  that  the  state  partial  differential  Eq.  (3-6)  and  the  costatc  partial-differential 
Eq.  (3-7)  may  l)o  soh'od  over  each  of  these  z  intervals,  for  each  oi,  as  a  constant  coefficient 
linear  differential  equation.  It  is  this  simplification  that  makes  computational  solution  of  the 
necessary  conditions  practical;  the  state  and  costate  may  be  propagated  from  the  beginning  to 
the  end  of  a  constant-u  interval  of  z  merely  by  multiplying  the  appropriate  transition  matrix 

i\7.  “  A  *  z 

(e  '  or  e  ').  No  numerical  integration  of  the  differential  equations  is  necessary. 

Wo  now  examine  the  behavior  of  the  Hamiltonian  over  a  constant-u  interval  of  z: 

H(x,p,u,z)  \  p'(z,  w )  A  (w,  u(z)|  x(z,  w )  do)  —  A  I  u  —  u  1^/2  (3-19) 

for  constant  u(z),  and 

=  \  (~p'(z,  w)  A^(w,  u)  x(z,  w)  -f  p'(z,  w)  A^(w,  u)  x(z,  w)l  dw  =  0  ,  (3-20) 

‘'wen 

The  Hamiltonian  is  therefore  constant  over  each  interval  where  u  is  constant.  Further,  the 
condition  for  switching  between  u^^  and  is 

IKX'h  p-!',  U|^,  z)  =  H(x‘!-,p’>,  Uj^^^.Z) 

The  Hamiltonian  is  thus  continuous  across  a  time  when  the  control  switches. 

These  two  conditions  are  sufficient  to  guarantee  that  the  Hamiltonian  is  constant  for  all 
zc  [-f,0|. 

We  may  similarly  differentiate  the  inner  product  of  state  and  costate  in  a  constant-u  inter¬ 
val  of  z: 


J}_ 

Oz 


wen 


p'  (z,  w  )  x(z,  w )  dw 


wen 


(— p'Ax  -f  p'Ax)  dw  =  0 


(3-21) 


.Since  p(z,  w)  and  x(z,  w)  are  continuous  in  z,  the  inner  product  shown  above  is  continuous 
across  a  location  where  the  control  u  switches. 

These  two  conditions  are  sufficient  to  guarantee  that  the  inner  product 


\  p'  (z,  w  )  x(z,  w )  dw 
*^wcn 

is  constant  for  all  zc  [—1,  0]. 

At  z  ■=  — /,  this  equals 


(3-22) 


‘wen 


K'(w)  x(— 1,  w)  dw 


(3-23) 
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but 


.I(x,  0)  i  K'(w)  X  (-<,  w)  da  +  ^  (  |u-u 
'■ujcV. 


dz 


Thus 


J(x  ,  0)  =  ^  p'(z,  w  )  x(z,  u>)daj+^r  |u— u  j^dz 

°  ''wen  "  ^ 


for  all  zf  (— /,  0],  in  particular  for  z  -  —0.  This  yields 


(3-24) 


(3-25) 


j(x  ,  0)=r  p'{0,  w)  X  (w)  dw  +  V 

°  ^’wen  °  “  '  -f 


(3-26) 


We  may  consider  the  initial  location  and  state  as  a  pair  of  independent  variables  and  sup¬ 
press  the  subscripts  in  the  equation  above.  Then 

(3-27) 

and 


3J( 


Kx.z)  ^  r 

'^wen 


—  p’(z,  w  )  A(w)  x(w )  dw  +  -^  I  u  -  u^  I  ^ 


This  equation  is  the  function-space  equivalent  of  the  Hamilton-Jacobi  equation: 

—  =  Hlx*,  p<‘,  u"!',  z)  =  max  II(x*,  j)*,  u,  z) 

ucU 


(3-28) 


(3-29) 


F.  SUMMARY 

In  the  systems  considered,  the  linear  partial  differential  equation  of  the  state  plus  the  con¬ 
straint  that  the  control  lies  in  a  discrete  set  permits  us  to  solve  a  linear  constant-coefficient  dif¬ 
ferential  equation  for  each  w  over  each  interval  of  z  where  u  takes  on  a  constant  value.  A 
similar  sot  of  equations  holds  for  the  costate,  so  we  may  propagate  the  state  and  costate  through 
any  constant-u  interval  of  z  by  a  simple  matrix  multiplication  with  the  transition  matrices. 

This  fact  makes  it  possible  to  eliminate  all  numerical  integration  of  the  state  and  costate 
during  iterative  solution  of  the  necessary  conditions  of  Eqs.  (3-6)  to  (3-9),  and  makes  this  nu¬ 
merical  solution  computationally  feasible. 

To  apply  these  results  to  the  specific  examples  of  the  transmission  line  and  the  waveguide, 
we  need  only  specify  the  A  matrix,  the  K  vector,  and  the  physical  parameter  corresponding  to 
u’!'  (z). 
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IV.  TRANSMISSION  LINE  (TEM)  PROBLEM 


A.  INTRODUCTION 

The  problem  of  interest  is  wide-band  impedanee  matching.  We  are  given  a  wide-band  source 
of  radio  frequency  energy  and  a  load  impedance,  as  shown  in  F'ig.  I-l. 

The  load  impedance  is 

Zj^(w)  =  +  jXJw)  (4-1) 

wliicli  is  a  complex  function  of  frequency  defined  over  the  frequency  band  of  interest,  uen. 

The  source  impedance  is 

Zg(w)  =  Rg(w)  +  jXg(a))  (4-2) 

(another  complex  function  of  frequency  for  wcR)  and  this  source  has  an  available  power  density 
S(a)),  a  real  function  of  frequency  also  defiticd  for  weR. 

Practical  examples  of  such  elements  might  bo  a  vacuum  tube  RF  amplifier  as  the  source  and 
an  antenna  as  the  load;  we  would  then  like  to  design  a  lossless  coupling  structure  or  impedance 
matching  structure  which  maximizes  the  real  power  delivered  from  the  source  to  the  load. 

The  first  example  of  such  an  impedanee  matohing  structure  is  the  transmission  lino  or  TEM 
structure  shown  in  Fig.  IV-1.  Mere  the  structure  length  is  fixed  as  /.  The  transmission  line  is 
considered  lossless  so  it  is  completely  charaoterlzed  by  its  inductance  per  unit  length  L(z)  and 
its  capaeity  per  unit  length  C(z).  Note  that  both  these  .are  dependent  on  position  z,  but  not  on  the 
frequency  w. 


LOSSLESS  TRANSMISSION  -IN'-  _ 

CHARACTERIZED  BY  L(2l,  C(2) 

Fig.  IV-1.  Impedance  match  using  o  TEM  structure. 


By  standard  convention,  real  power  flow  is  in  the  +z  direction,  so  we  connect  the  source  at 
z  =  — f  and  the  load  at  z  =  0.  Since  the  source  has  available  power  density  S(w),  we  consider  that 
a  forward  wave  w'ith  power  density  S(a,’)  flows  from  — f  to  0,  where  partial  reflection  occurs  and 
a  reflected  wave  having  power  density  T(a;)  returns  from  0  to  —  f.  The  net  power  delivered  to  the 
load  is  then 


V  = 


T(w)l  dw 


(4-3) 
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The  line  is  lossless,  so  S(w)  anil  T{a')  are  naturally  inilependent  of  z.  The  power  "lost" 
through  reflection  is 

.r  =  \  Tlca)  da'  .  (-I--!) 

In  the  design  procedure,  we  may  eithei-  maximize  \'  os-  minimize  J;  both  ai  e  exactly  equiv¬ 
alent  because  the  line  is  lossless.  The  problem  we  consider  is  thus  to  select  L(z)  and  C'(z). 

7.  |-f  .  0].  subject  to  appropriate  constraints  in  such  a  way  as  to  minimize  .1, 

The  particular  constraints  selected  for  this  design  are  the  following; 

(t)  L(z)  =  L  is  constant  for  all  z  (  t-'i) 

(2)  C(z)  f  (Cj,  C2 . ^m^  where  0  <  <  ■«  for  all  i  (-l-b) 

These  constraints  were  selected  to  make  the  physical  implementation  of  the  resulting  struc¬ 
ture  practical.  Consider  a  stripline  formed  by  parallel  conducting  strips;  to  synthesize  a  trans¬ 
mission  line  having  constant  inductance  per  unit  length  and  piecewise  constant  capacity  pei'  unit 
letigth.  one  need  only  fill  the  space  between  the  conductors  with  successive  intei'vals  of  vai'ious 

dielectric  materials  (see  Fig.  IV-2).  The  available  capacities  {Cj.  C ^ . should  be  chosen 

to  correspond  to  the  capacities  realized  by  standard  or  readily  available  dielectric  materials. 


I  •  -/  I  «o 


Fig.  IV-2.  Stripline  reolizotion  of  TEM  structure. 


The  construction  utilizes  dielectric  materials  rather  than  magnetic  materials  because  magnetic 
materials  having  constant  magnetic  susceptance  over  a  wide  range  of  frequencies  are  not  readily 
available. 


B.  THE  STATE  PARTIAL  DIFFERENTIAL  EQUATION 


A  TEM  (transverse  electro-magnetic)  structure  or  transmission  line  is  c  haracteidzed  by  an 
inductance  per  unit  length  l.(z)  and  a  capacity  per  unit  length  C(z)  if  we  assume  a  lossless  line. 
Line  losses  are  accounted  for  by  Ii(z).  the  resistance  per  unit  length  of  the  conductors,  and  (i(z). 
the  conductance  for  unit  length  of  the  material  between  the  conductors. 

3  0 

Voltages  (V)  and  currents  (T)  on  the  line  satisfy  the  following  differential  equations; 


ay; 

oz 


av 

til - QV 

at  ^ 


(-1-7) 


where  z  is  the  spatial  distance  down  the  line  and  t  is  the  time  variable. 
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Assuming  sinusoidal  stoady-statc  with 


from  one  end  of  tlic  interval  to  the  other  using  a  single  matrix  nuiltiplieation  by  the  transition 
matrix  for  eaeh  co; 


x(z  +  cl,  co)  = 


I  i-i  Ll  I  1  —  QOS  Oi 

'  2  ■  '  C  '  2 


C  ,  1  —  cos  O, 
L  '  2  ' 


-  s  in  O 


,  1  +  cos  Oi 
'  ■>  ' 


-.Ijc 


./I. 

n/2C 


sinO  0 


Itl 


sin  O 


cos  O 


x(z,  to) 


(4-11) 


where  (1  =  2u\'  l.C  cl,  and  d  is  the  length  of  the  constant-capacity  interval.  This  eliminates  any 
numerical  integration  of  the  partial  differential  equation  and  makes  the  computation  practical. 

The  boundary  conditions  on  the  state  are  as  follows:  at  z  =  0  a  load  ^  ^^1  Pon" 

neeted.  Remembering  the  definition  of  impedance  in  terms  of  voltage  and  current, 


x(0.  0.')  = 


IvTo.o))!*^ 

|Z,(a>)|^ 

|I(0,  a')|^ 

1 

\SS^^{0,x) 

v'8.Sj(0.  a') 

'Tzx^ix) 

\'2Rj  (w) 

\^Sj^(0,  a') 

s^n^ix) 

(4-12) 


This  provides  three  conditions  on  the  state.  The  fourth  condition  appears  at  z  -•  f: 

|V(-f,cu)^ 


Six)  = 


8Rg(w) 


[1,  |Z„(a))|^.  'Jlxjx).  \l\\Aw)] 


The  refleeted  power  is  also  evaluated  at  this  point; 


T(a')=  [1.  !Z<,(tu)|^.  \'^S^{x).  -^l\(x)] 


X' 


IK-J'.u))!' 


■fss^i-l,  x) 


\8S^i-l,x) 


!v(-f.  w)!' 

|l(-f.te)|^ 


V  8Sj(—  I ,  x) 


^^8S^^{-^,  x) 


(4-13) 


K'(ce)  x(— t ,  oj) 


(4-14) 
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C.  APPLICATION  OF  MAXIMUM  PRINCIPLE 
1.  Necessary  Conditions 

The  transmission  line  pi  ohlem  above  satisfies  the  restrirtions  of  the  development  described 
in  See.  Ill,  namely,  a  lineai'  system  with  linear  terminal  cost,  so  we  may  restate  the  necessary 
conditions  developed  in  Sec.  ill; 


i>x  '  {■/..  u)  _  r51K  (x.ji,  C',  ^ 
!)■/.  ~  (5p  '{z.u'Y 


A(a.' ,  C'l  X  (z.  0,' ) 


(4-15) 


itp  (z.Cc')  hll  >(x,  p,  C,  ■/)  ,,,  V 

ll(.\  ,  |5’'’ .  C  ,  z)  >  Il(x  ,  p'!" ,  C,  z)  for  ail  Cr  ( t  j.  (  2 . ^ 


\vh>’re 


and  the  Hamiltonian  is 


.1  =  1  T(a;)  da.'  +  (  4  ^z 

4'co  J-t  2 

=  ^  K'(w)  x(— f ,  w)  da'  ^  ^  (C'  —  ( 


iKx,  [),  C,  z)  =  ^  p'(z.  w)  A(a',  C')  x(z,  w)  dw  —  ^  (C  —  C 

''Wf  n  ^ 


(4-16) 

(4-17) 


(4-18) 


(4-19) 


Note  the  addition  (as  in  Sec.  Ill)  of  the  small  term  qiiadriitic  in  capacity  (control).  Here 
denotes  some  intermediate  capacity,  such  as  +  Cj^jj.^,)/2,  where  the  MAX  and  MIN  sub¬ 
scripts  denote  the  extremes  in  . ^m^ ' 

The  boundary  conditions  on  the  state  and  costate  arer 


x(0,  u)  =  x^(a’) 


(4-20) 


p(- 


■f,w)  =  [x(-f,  w),-fl  =  ~  ^  K'(w)  x(-f,a')  dca  =  K(a') 


''u}tU 


(4-21) 


2.  Properties  of  Solutions 

The  necessary  condition  of  Eq.  (4-15)  is  merely  a  restatement  of  the  system  equation 
Eq.  (4-10).  Equation  (4- 16)  yields 


^  p(z,  w)  =  -  xfZu 


0  -L  0 


0  C(z)  0 


C(z)  -L  0 


0  0 


p(z,  w)  =  —  A'(w,  C)  p(z,  u) 


(4-22) 
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Over  any  constant-C  interval  this  is  a  linear  constant  coefficient  differential  equation  and 
may  be  solved  using  the  transition  matrix: 


I  ^  ^  _eos  1 
'  '  2  ■  ’ 


£  /  1  ~  eos  O  1 
L  '  2  ' 


p(7  +  d.  w)  = 


L  ,  1  —  cos  O  . 
C’  '  2  ' 


it 


/  1  cos  8  V 
'  2  ' 


sinO 


sin  O  0 


cos  O  0 


p(z,  w) 


(4-23) 


where  O  =  Zu  si  LC  d.  and  d  is  the  length  of  the  constant-capacity  interval.  This  eliminates  any 
numerical  integration  of  the  costate. 

The  boundary  condition  on  the  costato  is; 


p{— t ,  w)  =  K(aj)  = 


|Z„(a.)| 


8Rc(a)) 


(4-24) 


^/2XJ,(w) 

sriRglw) 

Note  the  duality  between  this  condition  and  the  state  boundary  condition.  Eq.(4-12).  In  fact, 
this  duality  goes  much  further:  imagine  opening  the  transmission  line  at  any  point  ze  [— f,  0). 
Some  impedance  will  be  measured  at  that  point  looking  to  the  right  (  +  z  direction);  call  this  im¬ 
pedance  7.^{7.,u)  =  R^iz.u.')  +  jX^(z,  w).  Similarly,  one  sees  an  impedance  Z_(z,w)  looking  to  the 
left,  toward  the  source  (see  Fig.IV-3).  In  general 

r  |Z^(z,w)|^1 


x(z,  w)  = 


v^X^(z,w) 


2Sp(z,  w) 
R^(z,  u) 


(4-25) 


‘JZH^iz,  oj) 
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r  ' '  ~)  0  •  i 


R_(j.ul  +  J  ■  Z_(»,IU  1  )  •  R,(7,u!  ♦  ix^(7,  oil 

M ► 

^  )  Q  - ^ 


BREAK 

IN 

LINE 


■|vp' 

|z,(».a,lp  ■ 

1 

X (2,W  )  ■ 

|i|' 

^s, 

1 

8P  ii.u) 

4 

-  ! 
.  ! 

|zh,.)|* 

yixq.iu) 

yes. 

yzRjj.uiJ 

-V^R_(7,i.«) 

Fig.  IV-3.  Impedance  interpretation  of  state  and  costate, 
by  definition  of  the  impedance  in  terms  of  voltage  and  current.  Similarly, 


p(z.  w) 


i 

8R^(z,  w) 


|ZJz.a,-)| 

4zXJz.u) 


-nTzRJz.  u) 


(-1-26) 


Thus,  the  state  is  a  characterization  of  the  impedance  seen  looking  to  the  right  at  any  point 
and  the  costate  charact'^rizes  the  impedance  looking  to  the  left.  If  we  evaluate  Eqs.  (4-25)  and 
(4-26)  in  the  special  cases  z  =  0  and  z  =  —  f,  respectively,  we  obviously  obtain  the  previously 
stated  boundary  conditions  on  the  state  and  costate. 

it  is  interesting  to  interpret  the  maximization  of  the  Hamiltonian  in  this  light.  If  we  elim¬ 
inate  the  added  cost  term  which  was  quadratic  in  capacity,  the  question,  "Is  there  some  allowable 
capacity  change  one  can  make  to  increase  the  Hamiltonian?"  is  exactly  equivalent  to  the  question, 
"Is  there  some  allowable  capacity  change  one  can  make  (over  some  dz  interval)  so  that  the  paral¬ 
lel  combination  of  Zjz.gi)  and  this  capacity  is  closer  to  Z*(z.u))?"  (Here  the  superscript  aster¬ 
isk  means  complex  conjugate.)  Either  of  these  two  problems  results  in  the  same  equation  and 
the  same  capacity  change. 

The  Hamiltonian  evaluates  to: 

H(x,  p,  C,  z)  =  —  y  (C  -  y  vTuiC  (Xj(z,  cu)  P2IZ.  +  Xj(z,  w)  Pjiz,  cu)]  dw 

—  \  vTwL  [x-(z,  w)  p,(z,  ui)  +  x^(z,  w)  p-(z,  u:)]  do;  (4-27a) 

-wen 

where  the  numerical  subscripts  denote  the  different  components  of  state  or  costatc  vector. 
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More  simply, 


iKx,  p,  C,  z)  =  -  4  (C  -  C  f  M  ('  +  constant  t 
z  a  c 


(4-27b) 


Tlic  capacity  multiplier  is  then 


M  =\i: 
c 


2  \  tj  (x-p,  + 


(4-28) 


As  shown  in  Sec.  Ill,  tliis  specifics  tlie  capacity  which  maximizes  the  Hamiltonian  as 

C«=rc{C,.r, . C  }  closest  to  C  +M  ''A  (4-29) 

'  1’  2’  rn  a  c 

It  is  possilile  foi-  'A  to  lie  precisely  half  way  between  two  available  capacities.  If 


tliis  condition  holds  over  an  interval  of  z,  we  have  the  singular  condition,  where 
dM 


d^Al 
c  _  c 


dz 


=  0 


dz 


These  relations  yield 


-’=2'’2’  ° 


(4-30) 


(4-3tn) 


and 


Ai 

dz^ 


=  —Z'TzL  |l  C  a)^(x,p,  +  x,p,)  do)  +  C  \  aj^(x,p, 

L  ''ufU 


+  XjP^)  da^ 


=  0  (4-31b) 


or 


r  ^ 

C  =  -Un  <^3Pt  "  >=2P3> 

Ln  ‘^^(>'3P2  +  ^lP3’‘^" 

If  either  Kq.  (4-3ta)  fails  to  be  true,  or  the  C  of  Ilq.  (4-3tr)  equals  neither  of  the  two  candi¬ 
date  values,  we  arc  not  on  a  singular  extremal  and  the  ambiguous  condition  only  lasts  over  an 
interval  of  z  having  measure  zero.  If  both  equations  can  he  satisfied  for  one  of  the  candidates, 
we  are  on  a  singular  extremal  and  the  capacity  indicated  in  Kq.(4-3tc)  is  optimal. 

3.  Summary 

We  have  now  specified  a  differential  equation  for  the  state  and  costatc  for  each  frequency  u, 
a  practical  method  of  solution  (the  transition  matrices),  a  complete  set  of  boundary  conditions, 
and  a  procedure  for  uniquely  specifying  C(z)  as  a  function  of  x(z,ctj)  and  p(z,cij)  even  in  the  case 
of  the  singular  extremal.  This  coupled  set  of  equations  may  be  solved  iteratively  on  a  digital 
computer,  and  when  convergence  is  obtained,  the  functions  obtained  will  satisfy  all  our  necessary 
conditions  for  optimality. 


t  This  constant  is  independent  of  C. 
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If  we  liad  not  included  the  modification  to  the  cost  by  a  term  quadratic  in  C,  we  would  have 
been  forced  to  follow  a  completely  different  approach;  namely,  to  obtain  the  solution  of  the  prob¬ 
lem  when  C  is  constrained  to  an  interval  by  using  Eq.f'l-Slc)  to  determine  C  over  all  intervals 
of  z  where  =  0,  This  solution  will  provide  a  smoothly  varying  function  C(z),  as  a  function 

of  z,  requiring  numerical  integration  of  the  state  and  costate  differential  equations.  Such  a 
procedure  is  computationally  very  slow.  Moreover,  after  convergence  is  obtained  we  are  still 
faced  with  the  problem  of  finding  some  piecewise  constant  approximation  to  C(z)  to  permit  con¬ 
struction  of  the  structure. 

The  procedure  selected,  on  the  other  hand,  has  the  benefits  of  computational  efficiency  and 
a  naturally  piecewise  constant  solution  C(z). 
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V.  WAVEGUIDE  (TE)  PROBLEM 


A,  INTRODUCTION 

In  the  dovclopment  cliscus.si’cl  in  tiiis  section,  we  arc  again  concerned  with  the  problem  of 
designing  a  lossless  impedance  matching  structure.  Ilowi'ver.  the  frequencies  of  interest  in 
this  case  are  generally  in  the  microwave  region,  far  higher  than  the  radio  froqueneies  dealt  with 
in  the  transmission  line  coupler. 

At  these  very  high  frequencies,  waveguides  are  used  rather  than  transmission  lines  for 
power  transmission,  because  the  waveguides  provide  consitlcrably  lower  power  loss  (dissipation). 
The  source  of  mierow'ave  power  might  be  a  klystron,  magnetron,  or  even  a  receiving  antenna. 

The  load  in  this  case  might  then  be  a  transmitting  antenna  or  the  input  circuit  of  a  UlIF  receiver. 

In  any  case,  the  design  goal  is  a  lossless  impedance  matching  W'avcguidc  structure  which 

maximizes  the  power  deli\  crcd  to  the  load,  or  equivalently,  minimizes  the  reflected  power. 

The  structure  to  bo  considered  is  a  uniform  rectangular  waveguide  whoso  axis  extends  in 
the  +z  direction  from  z  =  — i  to  z  -  0.  In  accordance  with  the  convention  of  real  power  flow  in 

the  +z  direction,  the  source  will  be  at  z  =  —I  and  the  load  at  z  =  0. 

This  waveguide  will  be  filled  (uniformly  across  every  cross  section)  with  lossless  material 
having  dielectric  constant  e(z)  and  magnetic  permeability  (i(z)  (Fig.  V-1).  The  design  job  is  to 
select  these  functions  under  appropriate  constraints  so  as  to  maximize  the  power  delivered  to 
the  load. 


Fig.  V-1.  Waveguide  structure. 


The  constraints  selected  are  the  following: 

(1)  p(z)  =  for  all  zc  [— I.O). 

(2)  €(z)  is  in  ■  ■  ■  « 

(3)  The  waveguide  is  operated  solely  in  the  TEj  g  tnode  and  every 
frequency  cuen  (the  frequency  band  of  interest)  is  above  the  low 
frequency  cut-off  correspxjnding  to 

The  first  two  constraints  guarantee  that  a  solution  for  the  dielectric  constant  must  be  a 
piecewise  constant  function  of  z.  The  structure  may  then  be  physically  realized  by  filling  the 
waveguide  with  successive  intervals  of  various  dielectric  materials.  Of  course,  ^2’  •  •  ■ 
should  be  selected  to  correspond  to  readily  available  materials. 

The  third  constraint  is  merely  a  guarantee  that  it  is  possible  to  propagate  real  power  down 
the  waveguide  at  all  frequencies  of  interest.  If  a  frequency  cu  is  above  cut-off  in  the  p^, 
material,  it  will,  of  course,  be  safely  above  cut-off  for  all  the  other  materials  considered. 
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B.  WAVEGUIDE  EQUATIONS^^ 

'I'he  waveguide  structure  under  consideration  is  shown  in  Kig.  V-1.  'I'he  cavity  is  filled  with 
material  whose  susceptance  e  and  permeability  p  varies  with  position  z  down  the  waveguide, 
but  at  a  given  z  is  constant  over  the  cross  section 

{xe  [O.a]}  X  {ye  [0,hl} 

There  arc  no  free  charges  or  currents  in  the  material,  so  Maxwell's  equations  yield: 

VxK  =  -p-5^  ,  VxTr=:c-^  .  (5-1) 

If  we  consider  the  transverse  electric  field  (TK)  solution. 


(5-2) 

(5-3) 

(5-4) 

(5-5) 

(5-6) 

(5-7) 

(5-8) 


Tins  solution  is  constructed  as  the  supcrpiosition  of  two  equal  amplitude  uniform  plane  waves 
in  phase  propagating  in  different  direction.;,  as  shown  in  Kig.  V-2. 

If  both  waves  are  equal  in  amplitude,  are  in  phasi',  and  propagate  at  an  angle  O  with  rc'spect 
to  the  X-axis,  where  cosO  =  X  2a,  then  the  two  waves  interact  destructively  and  cancel  to  zero 
along  the  entire  lines  x  =  0  and  x  =  a.  We  merely  set  the  solution  to  zero  for  x  <  0  and  x  >  a  to 
obtain  the  solution  for  waves  in  the  waveguide.  (Interpret  the  waves  inside  as  multiple  internal 
reflections.)  Superposition  of  the  two  waves  gives: 
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IMPEDANCE  71 


■  E  t 


|(ul-fl») 


WAVENUMBER  $  E/»I 


Fig.V-2.  Superposition  of  plone  Fig.V-3.  Equivolent  transmission 

woves  in  woveguide.  line  model  for  waveguide. 

j/3  1.  sin  O 

f  = -T  [2jE  sin(/3^X  •  cos  0)1  e  ° 

-  2F  —  -  j/3  z  sin  O 

H  =  — ^  [1  cos  O  cos  X  -  cos  O)  -  i  j  sin  O  sin  (/3  X  ■  cos  0)1  e 

TJ_  Z  O  X  O 

o 


(5-10) 


where 

■ 

Ey  ^  gives  the  Poyntlng  (power  flow)  vector  component  in  the  x-directlon.  It  is  purely 
imaginary,  indicating  that  power  merely  " rattles  back  and  forth"  in  the  x-direction.  However, 
fcf  X  rf  gives  us  real  power  flow  in  the  z-direction  (along  the  waveguide)  so  we  will  focus  atten- 

y  * 

tion  on  these  two  components. 

Define  r)  =  "radiation  resistance"  =  E„/I1„: 

y  X 


_  ^o  _  1  _ _ p 

sinO  ”  ve  '  sin 6  “  / - ; 


(5-11) 


2  2  2 
pe  —  jr  /a  lu 


Define  P  =  2rr/wavelength  in  z-direction. 


P  =  P^  sinO  =  \( ep  slnO  =  «/ pc  —  rr^/a^cu^ 


(5-12) 


This  gives  rise  to  the  equivalent  transmission  line  model  of  Fig.  V-3,  where  Ky  is  the 
analog  of  line  voltage  and  11^  is  the  analog  of  line  current.  We  have  the  two  "transmission  line" 
equations. 


dE  dll  .  /  ,  X 

-d  =  -  ^y 


(5-13) 


The  electric  and  magnetic  fields  above  are  continuous  in  z  and  could  be  used  as  the  system 
state.  However,  it  is  far  more  convenient  to  select  the  following  state  vector  x(z,tu): 
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x(7.,a'  ] 


sfS  S,(z,a) ) 


^/8’  Sj^  (z.o: ) 


(5-14) 


where  S|(z,a))  is  the  imaginary  part  of  the  Poynting  vector  Sj(z,u.')  =  Im  {!•;  ,II^*/2}  and  Sj^  is  the 

real  part  S  ,(z,tj)  =  Ho  (K  i!  ’>/2}. 
it  y  X 

This  choice  of  a  state  vector  yields  the  following  state  partial  differential  equation  for 
lossless  waveguide. 


^x(z,w)  =  vTw 


0  0  -u  0 


0  0  e(z,a,’)  0 


e(z,ui)  0  0 


0  0 


x(z,u))  =  A(ti!,c)  •  x(z,Ci))  (5-15) 


where 


e(z,w)  =  e(z) - 


2  2 


^oa  w 


(5-16) 


Since  c(z)  is  a  piecewise  constant  function  of  z,  Eq.  (5-15)  is  a  constant  coefficient  linear 
differential  equation  for  each  w,  which  may  bo  solved  over  a  constant-c  interval  by  using  the 
transition  matrix: 


X(Z  +  d,  U!)  = 


1  *■  cosO. 


e  ,  1  —  cos  O, 


sin  0 


0,1—  cos  O 


-) 


1  +  co.s  O 


2e 


sinO 


-./^sinO  0 


sin  O  0 


cos  O 


.  .\(z,  u.')  (5-17) 
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where  e  is  defined  by  Kq.  (5-16), 

O  -  2wy^d 

and  d  is  the  length  of  the  constant-f  interval,  This  eliminates  any  numerical  integration  of  the 
partial  differential  (-quation. 

The  boundary  condition  at  z  •  0  is  determined  from  the  "wave  impedance" 


h,  (w) 


K  (z  =  0,0,' ) 


H^(z  0,w) 


=  Rj  (o.  )  +  jXj  (oj) 


(5-18) 


which  is  a  measured  characteristic  of  the  load. 

A  similar  wave  impedance  is  available  for  the  source; 

r)<,(o,’ )  =  Rg(a})  +  iXg(w  ) 


(5-19) 


and  this  source  has  available  power  density  S(o)). 
The  boundary  conditions  on  the  state  arc 


|Ky(0.W)|‘ 


Ul(w)| 


x(O.w)  = 


|llx(0,aj)| 


nTS  S|(0,w) 


S^(0,w) 


(5-20) 


'fS  S,^(0,ai) 


V  2  ) 


This  provides  throe  conditions  on  the  state.  The  fourth  condition  appears  at  z  =  -/ 

|Ey(-i.W)|^ 

=  8R^)ll.hg(w)|^.  ■/ZXgiw),  Nr2Rg(w)l  |Hx(-^w)| 

nTs  S,(-i,w) 

^  Sp(-/,w) 


(5-21) 


T(w),  the  reflected  power,  is  also  evaluated  at  this  point: 


|Ey(-i,w)|‘ 


=  8n^ll,hs(w)|^.  •sr2X,,(w),  -sr2R,-(w)l 


\f  &  Sj(-/,w) 


'/8’Sp(-/,w) 


K'(w)  X  i-l.u] 


(5-22) 
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C,  APPLICATION  OF  THE  MAXIMUM  PRINCIPLE 
1.  Necessary  Conditions 

The  neeessar'  contlitions  for  optimality  arc  those  developed  in  See.  Ill,  namely;  Given  the  cost 


.1  ^  (  r{u))  dw  +  ( 


dz 


=  j  K'(w)  xl-f.cj)  dw  4  A  ^  (c-c„)^dz 


^wen 


(here  c  denotes  som(;  intermediate  dicleetrie  constant),  we  define  the  Hamiltonian  as 

cl 

II(x,p,e,z)  =  (  p'(z,aj)  A(a!,c)  x(z,w)  dai  -  y  (c  -  e 

Note  the  addition  (ns  in  See.  HI)  of  the  small  quadratie  term  in  e. 

Then  there  must  exist  a  costate  p’''(z,w)  sueh  that; 

,,  ,  , 

‘  6p’'(z.w) 


and 


H(x^  p>^',e*,z)  ^  IKx*, P’1‘, c,z)  for  all  c  in  {£^,€2- 
1  ho  boundary  oondittons  on  the  state  and  costate  arc: 


x(0,w)  =  x^{w) 


(5-23) 


(5-24) 


(5-25) 

(5-26) 

(5-27) 


.  A  -  6J  [x(-i,w),-ii  _  6  r 

p(-<,w) - - ^  \ 


K'  (w  )  x(-l,w  )  do:  =  K'  (w  ) 


(5-28) 


2.  Properties  of  Solutions 

Equations  (5-23)  are  merely  a  restatement  of  the  system  equation,  Eq.  (5-). 5).  Equation 
(5-24)  yields: 


az 


=  p(z,w)  =  —'fZu 


e(z,w)  0 


e(z,w)  0  0 


.  p(z,w)  =  -A'(w,e)  ■  P(z,w)  (5-29) 


where  e(z,w)  is  given  by  h'q,  (5-16). 
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This  may  be  solved  over  any  constant-c  interval  by  the  transition  matrix: 


p(z  +  d.  oj)  = 


j  l_+_COS_0j 


-2  I  ^  ~  cos  9. 
e  '  2  ' 


^  sinG 
2e 


e  ,  1  —  eos  ©1 


1  +  eos  G 


sin© 


sin©  0 


sin© 
2e 


eos  © 


p(z,  w)  (5-30) 


where 


©  =  Zu^Tcu  d 
o 


and  d  is  the  length  of  the  eonstant-c  interval.  This  eliminates  numerical  integration  of  the 
costate. 

The  boundary  condition  on  the  costate  is 


p(-<,  w)  =  K(w)  = 


1 

8Rg(w) 


|ns(co)|^ 
»/T  Xg(w) 


-n/T  Rjj(co  ) 


(5-31) 


An  interpretation  of  the  state  and  eostate  analogous  to  that  of  Sec.  IV  is  possible  here.  The 
state  is  a  characterization  of  the  wave  impedance  looking  to  the  right  (toward  the  load)  at  any 
point  in  the  waveguide,  and  the  eostate  is  a  characterization  of  the  wave  impedance  loolcing  to 
the  left  at  that  point.  Equation  (5-31 )  is  an  obvious  special  case  of  this  for  z  =  -t. 

The  Hamiltonian  evaluates  to: 


H(x,p,e,z)  =  -  y  (e  -  c  )^  +  e  \  -IZu  (x,p,  +  x^p,)  dw 

C  rr  ,  "  ’'l P3> 

^^(XjPi  +  X2P3)  do; - - 


dw 


(5-32) 


The  dielectric  multiplier  in  the  Hamiltonian  is  then 
M  =  \  tj(x,p,  +  X,  p,)  dw 


(5-33) 
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As  shown  in  Sec.  Ill,  this  specifics  the  capacity  which  maximizes  the  Hamiltonian  as 

£>!■  =  c  in  {e,,  e,,  .  .  .  e  }  closest  to  c  -  M  'a  (5-34) 

1  ^  rn  Q  c 

It  is  possible  for  ^A  to  lie  precisely  half  way  between  two  available  dielectric 

constants.  If  this  condition  holds  over  an  interval  of  z,  we  have  the  singular  condition,  whore 

dM  d^M 

_ ^  _ _ ^  _  A 


These  relations  yield: 


'(iicn 


u!  (XjPj  -  ’^2^2^  do)  =  0 


(5-3S) 


^'o^^  LcU  ‘^<’^3^2  "  ’‘lP3>  LcQ  "  ’‘2P3> 

c  = - - -  - - _  .  (5-36) 

Len  “^<’‘3P2  +  ’'1P3) 

If  Kq.  (5-35)  fails  to  be  true,  or  the  e  of  liq.  (5-36)  equals  neither  of  the  two  candidatf*  values 
for  c,  wc  are  not  on  a  singular  extremal  and  the  ambiguous  condition  only  lasts  over  an  interval 
of  z  having  measure  zero.  If  both  equations  can  be  satisfied  by  one  of  the  candidates,  wo  are 
on  a  singular  extremal  and  the  dielectric  constant  indicated  by  Mq.  (5-36)  is  optimal. 

D,  SUMMARY 

A  differential  equation  has  boon  specified  for  the  state  and  oostato  for  each  w,  as  well  as  an 
efficient  method  of  solution  (the  transition  matrices),  a  complete  set  of  boundary  conditions,  and 
a  procedure  for  uniquely  specifying  e(z)  as  a  function  of  x(z,w)  and  p(z,u))  even  in  the  case  of  the 
singular  extremal.  The  resulting  set  of  equations  is  almost  identical  to  that  obtained  for  the 
transmission  lino  problem  with  the  obvious  exception  of  variable  names.  This  wavt:guide  problem 
may  thus  conveniently  be  solved  by  the  same  iterative  digital  computer  scheme  (program)  which 
solves  the  transmission  line  problem. 

When  convergence  is  obtained,  the  functions  obtained  will  satisfy  all  our  necessary  condi¬ 
tions  for  optimality. 

As  in  the  transmission  line  problem,  if  wc  had  not  modified  the  cost  functional  by  the  term 
quadratic  in  c,  we  would  have  been  forced  to  use  Kq.  (5-36)  over  all  intervals  of  z  where 
M  (z)  =  0  to  determine  c*(z).  This  would  require  the  (computationally  inefficient)  numerical 
integration  of  the  state  and  costatc  equations,  then  the  approximation  of  this  continuous  e(z)  by 
some  piecewise  constant  function  to  permit  construction  of  the  structure. 

On  the  other  hand,  the  present  method  is  computationally  efficient  and  provides  a  solution 
€(z)  which  is  naturally  a  piecewise  constant. 
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VI.  NUMERICAL  PROCEDURE 


A.  INTRODUCTION 

'I’hc  nc'ccssiiry  conditions,  liqs.  (3-6)  to  (3-8)  and  (3-l6)  wcii'  proftrammccl  for  iterative 
solution  on  a  digital  computer.  Hejuation  (3-16)  is  u.sed  to  find  the  control  u(7.)  (('(/.)  for  tiu'  trans 
mission  line,  and  c{7.)  for  the  waveguide)  in  terms  of  the  state  and  costate  vectors  x(7..  u)  and 
p(x.,  ui).  .Since  the  same'  pi’oc('dui-es  ai’c  eari'icd  out  in  the  transmission  line  case  and  the  wave¬ 
guide  case,  both  problems  art'  solved  l)y  tlie  same  progi-am. 

Flow  eliarts  for  this  solution  procedure  are  given  in  Appendi.x  15,  and  a  complete  i^’ortran- 
IV  soui  ee  listing  is  provided  in  Appc  ndi.x  Once  Kq.  (3-l6)  is  us(>d  to  eliminate  the  control 
from  the  state  and  costate  equations,  W('  luivc'  two  coupled  four-v<'ctor  systems  to  solve  subject 
to  split  boundary  value  conditions.  This  solution  proci'eds  iti'ratively  by  guessing  the  function 
.x(  — f,a!)  [p(  — f,  ui)  is  given  in  Fq.  (d-Zd)  or  F.tp  (5-31)]  and  inti.'grating  the  equations  in  thi'  +7.  di¬ 
rection  to  7.  r  0.  Tlu're  .x(0,aj)  will  not  meet  the  boundary  cotidition  imposed  by  tlu'  load,  so  the 
initial  guexss  for  .x(— f,a.')  must  be  changed  and  this  process  itc'i’ated  until  .x(0,  w)  does  meet  th(’ 
load  boundary  conditions.  \t  this  time,  tlu'  solution  x(7..ui)  ()(/.,  ic)  and  11(7.)  will  satisfy  all  our 
necessary  conditions. 

it  is  interesting  to  note  that  at  each  iteration  we  have  obtained  a  solution  to  some  optimal 
control  problem,  namely  that  of  matching  to  the  load  which  corresponds  to  our  inti'rmediate 
result  for  x(0,  tv  ). 


B.  SIMULTANEOUS  SOLUTION  OF  NECESSARY  CONDITIONS 

As  noted  above,  we  start  at  7.  =  -I  with  p(  — f,  w)  given  by  Kq.  (d-20)  or  Kq.  (5-31)  and  an 
initial  guess  for  x(- . ,  w).  Wo  then  use  Kq.(d-28)  to  dctormiiio  M(-f)  (the  midtiplior  in  the 
Hamiltonian).  This  will  specify  a  starting  value  for  u(  — <)  (Kij.  3-16),  depending  on  which  thresh 
old  values  M  lies  between.  By  differentiating  the  expression  for  M(7,)  with  respect  to  7.,  we 
obtain  two  frequency  integrals  of  products  of  state  and  (-ostate  components  which  specify 

and  —  .  (K'qs.  (d-31a)  and  (d-31b)  or  (5-35)  and  (5-36)1 

d/ 


Kvaluate  these  terms  at  z  =  —I,  fit  an  osculating  |)arabola  to  the  M(7.)  function  at  this  point  and 
predict  at  what  7.  location  M(7.)  will  first  cross  one  of  its  threshold  values  (see  Fig.  \'i-l ). 

Now  that  we  have  a  guess  at  the  loc  ation  of  the  first  switch  point,  use  the  transition  matrices 
to  |)ropagate  x(-f,w)  and  p(  — £,w)  down  to  this  predicted  switch  point.  This  is  accomiHishcd  by 
the  subroutine  "STKP"  flow  ehai’ted  in  F'ig.  B-5.  Here  we  re-evaluate 


M(7.) 


(IM 

d7. 


and 


If  M(7.)  has  not  ycd  reachi'd  the  threshold  value,  fit  another  parabola  and  repeat  the  proe('ss. 
If  M(7.)  has  passc'd  its  threshold,  we  have  overshot  the  switch  point  and  a  parabola  is  fitted  to 
take  us  back  to  the  switch  |5oint.  This  process  repeats  itself,  as  indicated  in  flow  charts  15-2 
and  15-3,  convciging  ciuadratically  fast  on  the  switch  point.  (Throe  iterations  are  almost  always 
sufficient  to  converge  on  the  switchpoint  with  10  ^  accuracy.) 
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M(l) 


Fig.  VI-1.  Extropolotion  by  fitting  osculoting  porobolo  to  M(z). 

Oico  this  switch  point  is  determined,  the  subroutine  "SWITCH"  is  called  (Fig.  B-4)  to  change 

2  2 

to  a  new  eonstant-capacity  ititerval.  In  this  routine,  the  new  value  of  u  is  computed,  d  M/dz 

is  recomputed  (since  it  is  a  function  of  u),  and  a  parabola  is  fitted  to  carry  us  forward  to  the 

ne.xt  switch  point.  This  process  continues  to  tlic  end  of  the  lino,  (Fig.H-2),  generally  taking 

40  or  50  steps  in  a  line  one  wavelength  long  and  having  five  admissible  control  values. 

A  number  of  other  considerations  are  taken  into  account  in  the  actual  program,  of  course; 

sometimes  it  is  impossible  to  fit  a  parabola  witich  intersects  a  threshold  value  in  the  correct 

direction.  If  this  happens,  a  linear  fit  is  tried  and  if  that  fails  a  step  equal  to  about  1/20*'^  of 

the  line  length  is  taken  in  the  appropriate  direction.  Since  the  parabolic  and  linear  fits  some- 

th 

times  give  ridiculously  large  distance  estimates,  in  no  case  is  a  step  larger  than  l/20'^"  of  the 
line  length  taken.  Also,  one  is  not  permitted  to  step  off  either  end  of  tlie  line. 

This  process  now  provides  a  fast  and  efficient  mapping  from  the  initial  guess  to 

the  state  x(0,  w).  This  mapping  will  scive  as  a  single  step  in  the  larger  Iteration  process  de¬ 
scribed  below. 

C.  TCHEBYCHEFF  APPROXIMATIONS  FOR  BOUNDARY  FUNCTIONS 

In  the  iH'oblcms  considered,  the  state  and  costate  have  been  expressed  as  the  continuous 
functions  of  w  and  the  llai.^iltonian  was  expressed  as  an  integral  over  w.  In  actual  practice, 
howc-ver,  a  digital  computer  cai.not  handle  an  entire  function  numerically.  At  best,  one  can 
establish  a  grid  on  ui  and  approximate  the  integral  over  a;  as  a  sum  over  the  grid.  This  was 
done  in  the  program,  selecting  11  disci  ste  values  for  ai .  The  state  and  costate  transition 
matrices  wore  used  only  with  these  11  values  of  w,  and  all  ui-integrals  became  sums  over  these 
1 1  points. 

.Still,  roprf'sentation  of  the  source  impedan.''c,  load  impedance,  or  impedance  looking  into 
the  line  [whiidi  determines  our  x(-f,u))  guess]  cacli  required  22  numbers;  the  real  part  and  the 
imaginary  part  of  the  impedance  over  11  frequencies.  The  load  impedance,  the  guess  on  the 
initial  impedance,  or  the  resulting  impedance  at  the  other  end  of  the  lino  can  thus  be  considered 
22-vectors.  Iteration  to  convergence  in  this  22-spaec  seemed  numerically  unattractive,  so  it 
was  decided  to  use  a  more  compact  (lower  dimensional)  representation  of  these  impedance 
functions  -  the  Tchebyctieff  approximation. 
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ICtich  impedance  function  (defined  over  11  points)  may  be  considered  as  a  sampled  version 
of  a  linear  combination  of  the  first  11  Tcbcbychoff  polynomials.  These  polyn./mials  wore  se¬ 
lected  l)ecuuso  an  approximation  of  this  type  tends  to  minimize  the  maximum  approximation 
error  over  an  interval  (well).  I'urther,  a  series  expansion  of  a  function  in  Teliebycheff  poly¬ 
nomials  tends  to  converge  far  faster  (in  that  interval)  than  .some  otlii'r  polynomial  expansion;  it 
can  be  truncated  sooner  for  tlie  same  error. 

In  tlie  pr-olilcm  at  li.ind,  tlie  Tchebyeheff  expansion  was  truncated  after  4  terms,  resulting 
in  a  cubic  approximation  for  the  I'cal  part  and  another  for  the  imaginary  part  of  eacli  in.pedance 
function.  .Sueli  a  function  is  now  described  by  8  numbers  insti>ad  of  22,  and  iteration  proceeds 
in  an  8-dimensional  space.  'I'lie  Tclieljyclieff  functions  used  for  this  approximation  are  shown 
in  I'ig.  \T-2. 


Fig.  VI-2.  First  four  Tchebyeheff  polynomials. 


'I  h-<i-tiitoi 

To(.)  •  I 

0 - I 


The  real  part  of  a  load  impedance  is  repri'Si’Uted  by  11  numbers  (its  valui;  at  11  frequencies). 
To  convert  to  four  numljci'S  (tlie  first  4  coefficients  of  its  ’|■etlebyclleff  series  expansion)  it  is 
necessary  to  multiply  tlie  original  11 -vector  liy  the  4x11  nKitrix,  C^UAD,  described  below. 
l''ir.st,  form  tlie  11  x  11  matrix  A  sucti  that 

A..  =  T  ((i-6)  5] 

where  re|)r('sents  tlie  (J --  1)*’'^  Tclieliycheff  polynomial. 
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TC'IIKH  =  A 


(;in  1 1  X  •!  nintrix,  actually  a  truncated 
version  of  A) 


QL'AD  = 


(a  -1  X  H  matrix  which  is  a 
truncated  version  of  A"^) 


A  Tchebychoff  (luadrature  is  performed  by  prcunultiplying  an  11-vector  by  QUAD,  obtaining 
the  -l-vector  Tchebycheff  representation.  The  opposite  conversion  is  obtained  by  premultiply¬ 
ing  the  -l-vector  Tchebycheff  representation  by  TCHKH  to  obtain  the  equivalent  11-vector 
(function  of  frequency). 

It  is  simple  to  verify  that 


QUAD  X  TCHKH  = 


,  the  4X4  Identity  matrix  , 


so  tlie  desired  Tchebycheff  quadrature  (with  truncation  to  four  terms)  has  indeed  been 
aecomplisliod. 

The  program  in  Appendix  C  accepts  these  two  matrices  as  input  data  before  the  beginning 
of  computation.  If  for  some  reason  one  desired  to  expand  the  load  impedances  (functions  of 
frequency)  in  a  different  polynomial  series  representation,  say  P^(x),  it  is  only  neeessary  to 


Ai-Pj-i[(i-6)/5)  , 

in  the  above  procedure  and  input  the  two  resulting  matrices  to  the  program  in  place  of  QUAD 
and  TCHEB. 


D.  ITERATION  ON  BOUNDARY  CONDITIONS 

Now  we  have  a  compact  8-vector  representation  for  impedanee  functions  and  a  programmed 
procedure  to  map  between  such  functions.  Starting  with  an  initial  guess  on  the  impedance  pre¬ 
sented  by  the  coupler,  map  this  into  an  impedance  at  the  other  end  and  express  as  an  8-vector. 
Repeat  the  process  for  eight  small  Tehebycheff  variations  on  the  initial  guess  (the  first  four 
Tehebyehefi  polynomials  as  variations  to  the  real  part,  then  the  last  four  as  variations  to  the 
imaginary  part  of  our  guess  impedance).  Form  the  8X8  sensitivity  matrix  describing  the 


variation  in  the  terminal  8-vector  with  respect  to  our  8  variations: 

aTBUFF. 

SNSMAT{i,  j)  =  9tgUESS~ 

J 

where  TGUESS  is  the  8-vector  Tchcjbyeheff  representation  of  our  starting  impedance  guess,  and 
TBUFF  the  representation  of  the  resulting  load  impedance  obtained  by  the  mapping  procedure 
of  Fig.  B-2. 

To  obtain  a  gradient  vector  in  this  8-spacc  describing  the  changes  in  TGUESS  necessary  to 
drive  TBUFF  to  TLOAD,  the  representation  for  the  actual  load  impedance  we  are  attempting  to 
match,  a  linear  mapping  is  assumed  between  TGUESS  and  TBLT'F;  tlien 

GRAD  =  SNSMAT"^  (TLOAD  -  TBUFF)  , 
the  desired  gradient  vector. 

If  the  mapping  of  Fig.  B-2  wore  exactly  linear,  changing  TGUESS  by  GRAD  to  obtain  tlie 
new  guess  would  drive  TBUFF  to  TLOAD  exactly  in  one  step,  thus  providing  imm('diate  con¬ 
vergence  of  our  boundary-function  iteration  scheme.  In  reality  the  mapping  is  not  linear,  so 
the  process  described  above  must  repeat  at  this  point,  iterating  until  convergence  is  obtained. 
This  procedure  is  flow  charted  in  Fig.  B-1.  This  Newton-Raphson  procedure  will  converge  quad- 
ratically  fast  when  the  mapping  from  guess  impedance  to  resulting  load  impedance  is  approxi¬ 
mately  linear  over  a  region  containing  both  the  starting  guess  and  the  final  (solution)  guess. 

E.  SUMMARY 

We  have  described  an  iterative  procedure  that  iterat('s  on  boundary  conditions  of  the  problem 
and,  when  it  converges,  provides  quadratically  fast  convergence.  Of  course,  this  is  only  a 
general  description  which  skips  over  the  myriad  of  tiny  details  and  special  cases  with  which  the 
programmer  must  concern  himself. 

For  the  sake  of  completeness,  the  entire  program,  written  in  Fortran-lV,  is  provided  as 
Appendix  C.  This,  plus  the  above  general  description,  should  permit  mechanization  of  the  same 
procedure  on  another  system  or  in  another  language,  such  as  ALX30L  or  PL-1. 
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VII.  LIMITATIONS  OF  THE  PROCEDURE 


A.  SOLUTION  FOR  AN  INTERVAL  CONSTRAINT  ON  THE  CONTROL 
1.  Introduction 

In  all  of  tho  problems  considered  to  this  point,  the  control  ufz)  (either  C(z)  or  f(7.))  has  been 
constrained  to  lie  in  some  discrete  set; 

uf  {u, ,  u,,  .  .  .  u  } 

‘  1'  2’  nr 

This  constraint  was  imposed  for  two  reasons:  first,  the  control  function  is  then  piecewise 
constant  in  z,  [lermitting  the  use  of  transition  matrices  to  propagate  the  state  and  costate  down 
the  line;  this  eliminates  the  necessity  for  numerical  integration  of  the  state  and  costate  partial 
differential  equations,  yielding  far  better  accuracy  in  considerably  loss  computer  time.  If  the 
control  v/ere  meriTy  constrained  to  lie  in  an  interval  U,  the  piecewise  constant  property  of  the 
control  function  would  be  lost  and  numerical  integration  would  l)ecome  a  necessity. 

The  second  reason  for  the  difj^n'ote  set  constraint  is  convenience  of  physical  realization. 

An  interval  constraint  would  yield  a  smoothly  varying  f(z)-,  however,  a  material  having  this  same 
f(z)  is  not  readily  obtainable  and  would  be  extremely  difficult  to  synthesize.  One  would  wind  up 
making  a  piecewise  constant  approximation  to  this  f(z)  by  sanawiching  very  thin  sheets  of  the 
various  dielectrics;  it  would  be  far  simpler  in  the  first  place  to  solve  for  the  best  piecewise 
constant  f(z)  by  using  the  discrete  set  constraint. 

However,  the  solution  to  the  problem  when  an  interval  constraint  is  imposed  is  quite  im¬ 
portant  as  a  basis  of  comparison  for  the  performance  obtained  when  the  control  is  constrained 
to  a  (more  restrictive)  discrete  set  having  the  same  maximum  and  minimum  value.  For  this 
reason,  we  will  examine  the  interval  constraint  problem. 


2.  Sufficient  Conditions 

We  consider  the  following  system: 

(a)  A  distributed  parameter  system  whose  state  obeys  the  partial 
differential  equation. 

"  A(a',u)  •  x(z,  oi)  ;  well 

x(0,  u)  =  x^(w) 

where  A(w,  u)  is  linear  in  u. 

(b)  There  is  a  cost  functional  on  the  control  u(z)  of  the  following  form: 

n  po  (u  —  u  ) 

J{x  (a),u(-/,  0)}  =  \  K'(w)  x(-/,  w)  du;  +  A  \  - dz 

"'-I 


(7-1) 

(7-2) 


(7-3) 


Wc  want  now  to  find  the  optimal  control  u(z),  ze  (—1,  0),  which  minimizes  the  value  of  the 
cost  functional  J.  This  minimization  is  subject  to  the  constraint  that  u(z)  is  in  the  interval  U 
for  all  Z€  (— /,  0). 

Define  the  Hamiltonian 


H(x,  p,  u,  z) 


J 


p'(z,  ai)  A(a',  u)  x(z,  a)  da  -  y  (u  -  u^)^ 


(7-4) 


Tho  sufficient  conditions  for  optimality  derived  in  Sec.  11  yield: 


,  A(z,u)x(z,c.) 


-Sp 


x(0,  a,')  = 


,  _A.(z,u)p/z,a.) 


p(-t.u)  =  K(u.') 

H(x’!' ,  p*  ,  u*  ,  z)  >  II(x*  ,  p*  ,  u,  z)  for  all  ut U 

with  equality  holding  only  for  u  =  u*  . 

The  h’rechct  derivative  of  the  optimal  cost  function  J(x,  z)  must  exist: 

6.1(x.  z) 

6x 


(7-5) 

(7-6) 

(7-7) 

(7-8) 

(7-9) 


is  a  valid  I'rechet  derivative,  i.e., 

lim  .llxfg.’)  t  h(w),  z)  -  J(x(a)).  z] 


/^^,j.nlh(a.')|dw-0 


■Ld} 


(7-10) 


exists  for  all  h(a.')  having  compact  support  in  fl. 

First,  examine  the  validity  of  Kq.(7-10)  in  the  present  system;  this  required  Frechet  de¬ 
rivative  certainly  exists  for  z  =  —I,  because  we  can  differentiate  Kq.  (7-3)  with  respect  to 
x(-f,  w)  to  obtain: 


6.T/6x  obeys  the  same  partial  differential  equation  as  p(— /,  w)  if  it  exists: 
Oz  ^31  -A'(w.  u)  [x(w),  z]) 


(7-11) 


(7-12) 


so  6J/6x  merely  obeys  a  linear  partial  differential  equation. 

Since  6.1 /fix  is  well  defined  at  the  boundary  z  =  —I  (and  is  there  equal  to  K(ai)],  it  will  be 
well  defined  for  all  zc  (— /,  0).  The  system  does  satisfy  the  sufficient  condition  Eq.(7-10). 

Now  check  the  validity  of  sufficient  condition,  F.q.  (7-9);  we  know  that  u(z)  is  independent 
of  w  and  A(w,u)  is  linear  in  u.  The  u  may  then  be  taken  outside  the  frequency  integral  in 
Eq.(7-4),  demonstrating  that  the  Hamiltonian  is  a  function  of  u  of  the  form 


11(  X,  p,  u,  z)  =  —  ^  (u  —  u^)^  +  M  •  u  +  constant 
Since  A  >  0,  there  is  a  uniquely  maximizing  value  of  ucU: 


(7-13) 


u*  = 


+  m/a 

a 

if  u  +  m/a  is  in  interval  U 
a 

‘MIN 

if  u^  +  m/a  <  every  ueU 

‘MAX 

if  u^  +  m/a  >  every  ueU 

(7-14) 
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Thus,  the  system  under  considerntinn  satisfies  sufficiency  condition  l’;q.(7-9).  The  re¬ 
maining  conditions  are  exactly  the  same  as  the  necessary  conditions  for  the  same  system;  if 
any  solution  u(z)  is  found  which  satisfies  the  (iK'cessary)  conditions  i:q.(7-5)  to  (7-9),  then  this 
solution  is  the  unique  globally  optimal  solution  under  the  constraint  ur  interval  U. 

Merely  solving  the  necessary  conditions  for  ucU  yields  the  unique  optimal  control  which 
will  be  the  basis  for  our  performance  comparisons.  Since  the  set  {uj,  u^,  .  .  .  ^  subset 

of  U,  the  value  of  the  cost  functional  for  the  discrete  set  constraint  must  be  ^  the  cost  func¬ 
tional  for  the  less  restrictive  ucT  constraint. 

In  the  limit,  as  the  set  {u^,  U2,  .  .  .  u^^}  is  allowed  enough  members  to  become  dense  on  the 
interval  U,  the  optimal  cost  functional  for  this  constraint  approaches  the  cost  for  the  interval 
constraint.  In  addition,  the  control  function  u*(z)  converges  in  the  1.^  norm  to  the  control  func¬ 
tion  for  the  interval  constraint  because  of  the  uniqueness  property  of  this  solution  and  the  smooth¬ 
ness  (Lipschitz  condition)  of  the  system  partial  differential  equations. 


B.  SOLUTIONS  OBTAINED  FOR  CONTROL  CONSTRAINED  TO  A  DISCRETE  SET 

1.  A  First  Variation  Test 

Let  us  suppose  that  by  numerical  means  we  have  found  a  solution  u(z),  x(z,  w),  p(z,  w)  for 
zf  [-<,  0|  and  wen  which  satisfies  all  the  necessary  conditions  Eqs.  ( 3-6),  (3-7),  (3-8),  (3-16) 
under  the  constraint  ue  {u^,  u^,  .  .  .  u^}  . 

This  control  u(z)  is  necessarily  a  piecewise  constant,  and  is  a  candidate  for  the  optimal 
control.  It  is,  of  course,  not  guaranteed  to  be  the  optimal  control  because  it  satisfies  only 
necessary,  not  sufficient,  conditions.  In  this  case,  the  sufficient  condition  of  Eq.(7-10)  (cost 
smooth  in  x)  is  satisfied,  but  the  requirement  for  a  u  which  uniquely  maximizes  the  Hamiltonian 
is  not  satisfied  whenever  two  u  values  are  equally  close  to  u^  +  M/A. 

First,  let  us  suppose  that  this  candidate  solution  does  not  correspond  to  a  singular  extremal, 
i.e.,  the 

dM  _  d^M  _  Q 


relations  were  not  explicitly  used  during  the  solution  of  necessary  conditions  of  Eqs.  (3-6)  to 
(3-8)  and  (3-16).  Now  we  can  interpret  the  results  of  the  maximum  principle.  Note  that  the 
control  variation  analyzed  to  derive  the  necessary  conditions  of  Sec.  11  is  a  'first  variation,  in 
particular,  a  narrow  pulse  added  to  the  candidate  control  function.  The  result  is;  if  the  nec¬ 
essary  conditions  of  Eqs.  (3-6)  to  (3-9)  are  satisfied,  such  a  variation  of  the  control  can  only 
increase  the  cost  J.  Two  examples  of  this  type  of  perturbation  on  the  control  are  shown  in 
F  ig.  Vll-1  applied  to  a  candidate  control  which  is  a  switching  solution. 

As  shown  in  Fig.  Vll-1,  the  specified  first  variations  may  merely  correspond  to  a  first 
order  (c-sized)  change  in  switching  time  t.  The  analysis  yielded  the  result  that  6J  >  0  for  each 
variation,  6J  =  0  for  e  =  0,  and  J  is  continuous  in  e.  These  combine  to  yield; 


M 

6t 


T* 


=  0 


and 


JlT*  )  ■$  J(t)  for  T  w  T* 


where  we  are  considering  the  switching  time  t  to  be  the  independent  variable  we  are  perturbing. 
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CANDIDATE  SOLUTION 


Fig.  VII-1.  Two  possible  first  voriotions. 

In  short,  if  a  candidate  switching  solution  u(t)  satisfies  the  necessary  conditions  of  Kqs.  (3- 
to  (3-9)  then  the  cost  J  is  at  a  local  minimum  with  respect  to  each  of  the  switching  times  Tj; 
any  perturbation  of  these  limes  can  only  increase  the  cost  .1. 

2.  A  Second  Variation  Test 

Now  we  consider  a  second  variation  of  the  control  u(t),  in  particular: 

0  for  t/  (t  -  c,  T  +  c) 
u(t)  =  a  for  tc  (t  -  c,  t) 

-a  for  tf  (t,  t  +  c) 

Such  a  perturbation  may  be  called  a  "doublet"  pulse. 

To  utilize  the  results  of  previous  work  dealing  with  second  variation  tests,  we  must  first 
define  a  modified  Hamiltonian: 

A 

ll(  X,  p,  u,  t)  =  convex  cover  of  ll(x,  p,  u,  t) 

The  function  H(x,  p,  u,  t)  is  a  linear  function  of  u  and 


over  intervals  of  u  where  H  is  undefined  (u/  {u^,  U2, . . .  •’forming  a  linear  interpolation 

between  points  u, ,  u,,  .  . .  u„  where  H  is  defined. 

12m  — 

Assuming  that  the  candidate  switching  solution  is  not  a  singular  extremal,  then  it  is  easily 

A 

seen  that  minimizing  H(x,  p,  u,  t)  with  respect  to  u  yields  precisely  the  same  switching  solution 
as  in  Kq.(3-16)  which  minimizes  IKx,  p,  u,  t).  This  requirement  may  be  substituted  for  Kq.(3-16) 
yielding  the  relation 


over  every  interval  of  u  where  switching  occurs  between  the  endpoints. 

31 

Now  we  may  use  the  results  of  H.  J.  Kelley  in  determining  the  second  variation  of  the  cost 
under  the  "doublet"  control  variation  described  above.  The  result  is 


.2,  223 

6  J  =  —  j  a  € 


A.  (jL 

au  \  .,2  du  / 


+  higher  order  terms  in  e 


(7-15) 


.as  long  as  the  perturbed  control  remains  in  an  interval  of  u  where 


Thus,  if 


=  0  , 


3  d*-  3H 

3S  n  <3n> 


no  second  variation  on  the  control  will  reduce  the  cost.  If,  however,  this  term  is  positive,  an 
allowable  second  variation  will  reduce  the  cost;  the  candidate  solution  is  not  optimal. 

An  allowable  second  variation  is  shown  in  Fig.  VTI-2;  note  that  the  effect  of  this  variation 
is  10  cause  a  "chattering"  transition  rather  than  a  single  switch.  The  cost  can  be  reduced  still 
more  (assuming  that 


at  each  of  the  three  switching  times  of  Fig.  VII-2)  by  making  additional  second  variations,  i.o., 
turning  each  of  the  three  transitions  into  a  3-transition  chatter  (see  Fig.  VII -3).  Similarly,  wo 
reduce  the  cost  by  going  from  9  to  27  transitions,  etc.  The  limit  of  this  process  is  not  well 
defined,  as  it  corresponds  to  infinitely  fast  chattering  over  some  interval  of  time,  in  particular, 
the  interval  of  time  where 


3  d*-  ,3H. 

35  I  as'l  ^  » 


However,  an  interpretation  of  this  control  behavior  is  still  possible  in  terms  of  duty  cycle: 
if  u(t)is  chattering  infinitely  fast  between  u^  and  U2,  taking  on  the  value  of  u^  a  fraction  K  of  the 
time  and  U2  a  fraction  (1  —  K)  of  the  time  (of  course,  these  fractions  must  be  defined  as  Lebesque 
integrals  on  the  measure-space  of  a  distance),  then  the  trajectories  of  the  state  and  costate  are 
identical  to  those  observed  when  the  control  takes  on  the  value  Kuj  +  (1  —  K)  U2  (Ref.  32). 
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I 


I 


Fig.  VII-2.  A  possible  second  vorioHon 


at  a  switching  point,  implies  that  the  eanclidate  solution,  whieh  makes  a  single  transition  at  the 
switch  point,  is  inferior  to  a  eontrol  whieh  takes  on  intermediate  values  around  the  switch  point  — 
a  continuous  change  from  Uj  to  U2  will  result  in  a  lower  eost.  This  smooth  transition  takes 
place  over  the  interval  where  Kq.(7-16)  holds. 

This  simple  test  [Kq.  (7-16)1  permits  us  to  evaluate  the  performanee  of  the  eandidate  solu¬ 
tion  by  eomparison  with  the  solution  obtained  in  See.  VllA  under  an  interval  eonstraint.  Let 
the  optimal  cost  for  this  interv  al  constraint  be  denoted  by  Jj,  and  the  (sometimes  eontinuously 
v’arying)  control  function  by  Uj(z). 

At  each  switching  point  of  the  candidate  solution,  if  relation  (7-16)  fails  to  hold,  then  the 
"doublet"  control  perturbation  at  that  point  will  merely  increase  the  eost  functional.  This  happenr 
over  those  regions  of  z  where  "bang-bang"  solutions  u(z)  appear,  switching  between  and 

^MAX'  says  that  around  the  switch  point  in  question,  u(z)  provides  an  optimal  approxima¬ 

tion  to  Uj(z)  under  the  constraints  imposed.  On  the  other  hand,  at  switch. points  where  Eq.  (7-16) 
holds,  there  should  be  a  continuous  transition  rather  than  a  step  discontinuity;  a  control  per¬ 
turbation  that  approaches  this  smooth  Uj(z)  transition  will  decrease  the  eost  because  <  0 
(see  Fig.  VH-4). 

Instead  of  approaching  the  Uj(z)  behavior  by  incorporating  infinitely  fast  "chattering"  (which 
is  impractical  for  a  digital  computer  numerical  scheme),  we  need  merely  relax  the  control 
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Cll) 


constraint  uc  {u^,  U2, . . .  admissible  control  values  •  •  •  “n^’ 

fall  bet'veen  the  previous  values.  The  cost  for  those  new  control  values  is  such  that 

H(x,  p,  -  2H(x,  p,  P*  ’^k'  ^  °  (7-17) 

where 

“k  <  “k+1/2  <  V+1  • 

A 

This  is  a  type  of  convexity  defined  over  discrete  points  of  u.  Again  talto  H(x,  p,  u,  t)  as  the  con¬ 
vex  cover  of  H(x,  p,  u,  t). 

Tile  algorithm  previously  used  will  now  converge  to  a  switching  solution  which  takes  on  the 
newly  defined  intermediate  values  of  control  over  intervals  of  time  where  before  there  was  a 
single  transition  between  Uj^  and  This  stairstep  behavior  is  shown  in  Fig.  VI-1,  and  this 

"interpolating"  behavior  is  a  consequence  of  the  fact  tnat  M(z)  is  a  continuous  function  of  z. 

The  total  cost  of  this  solution  will  be  less  than  the  previous  solution  because  a  second  var¬ 
iation  decreases  the  cost  as  long  as  Eq.  (7-16)  holds.  The  convexity  assured  by  Eq.  (7»17)  fur¬ 
ther  decreases  the  cost,  making  it  less  expensive  to  take  on  the  value 

“k+1  +  “k 
2 

than  chattering  infinitely  fast  with  50  percent  duty  cycle  between  Uj^^j  and  Uj|,.  However,  as¬ 
suming  that  Eq.  (7-16)  still  holds  at  each  of  the  two  switch  points  in  Fig.  VI-1,  the  cost  could  be 
reduced  still  more  by  providing  a  smooth  transition  rather  than  a  switching  for  each. 

The  same  argument  again  applies,  so  defining  a  still  finer  grid  of  jaermissible  u  values  and 
repeating  the  solution  of  Eqs.(3-6)  to  (3-8)  and  (3-16)  will  result  in  a  lower  cost. 

Making  the  partition  on  available  values  of  u  increasingly  finer  can  only  decrease  the  value 
of  the  cost  functional  J,  because  the  new  solution  corresjxinds  to  a  less  restrictive  set  of  con¬ 
straints  on  the  problem.  The  cost  obtainable  as  the  partition  becomes  finer  is  bounded  from 
below  by  the  cost  Jj  corresponding  to  the  solution  with  an  interval  constraint  on  the  control. 

Since  J  is  monotone  nonincreasing  and  bounded  from  below  as  the  fineness  of  partition  is 
increased,  it  musfapproach  a  limit  J* .  In  fact,  J*  =  Jj  and  u(z)  converges  (in  norm)  to 
vl(z),  the  unique  control  solution  for  an  interval  constraint,  as  the  points  of  available  capacity 


I  ecome  more  dense  in  the  interval  U.  (This  is  a  consequence  of  the  smoothness  (Lipschitz  con¬ 
dition)  of  the  system  partial  differential  equations,] 

This  result  provides  a  straightforward  means  of  evaluating  just  how  close  to  optimal  (in 
terms  of  cost)  a  particular  suboptimal  solution  is.  Having  obtained  the  solution  and  the  cost 
when  we  are  constrained  to  n  values  of  capacity,  repeat  the  solution  with  finer  partition  to  ob¬ 
tain  etc.  When  one  is  satisfied  that  the  cost  has  converged,  i.e.,  a  finer  partitioning 

of  available  capacity  values  will  not  result  in  a  cost  which  is  measurably  lower,  then  that  cost 
equals  J*  =  Jj.  The  cost  improvement  that  has  been  obtained  through  finer  partitioning  is  pre¬ 
cisely  the  amount  by  which  the  original  candidate  was  suboptimal. 

This  process  of  obtaining  solutions  for  successively  finer  partitions  on  u  may  be  viewed  as 
a  method  of  numerical  integration  of  the  differential  equations  which  correspond  to  solution  of 
the  necessary  conditions  with  an  interval  constraint  on  the  control.  In  this  case,  discretization 
occurs  for  the  interval  U,  rather  than  the  more  common  discretization  of  [— i,  0).  "  with  any 

scheme  for  numerical  integration,  the  burden  lies  on  the  user  to  decide  when  partitioning  is 
sufficiently  fine,  i.e.,  when  making  the  partition  still  finer  will  not  result  in  a  significant  change 
in  cost. 


3.  Discussion 

We  have  noted  previously  that  a  candidate  ^witching  solution  which  satisfies  I!qs.(3-6)  to 
(3-8)  and  (3-16),  minimizes  the  cost  with  respect  to  each  of  the  switching  times  any  per¬ 
turbations  of  these  times  can  only  increase  the  cost.  If  Hq.  (7-16)  does  not  hold  at  these  switch 
times,  even  a  second  variation  cannot  decrease  the  cost.  However,  at  those  switch  points  where 
Eq.(7-16)  holds,  the  cost  may  be  reduced  by  either  the  chattering  behavior  shown  in  Fig.  Vll-2, 
or  the  stair-step  control  behavior  shown  in  Fig.  Vll-3. 

Either  of  these  changes  increases  the  number  of  switchings  in  the  control  function.  If  there 
is  some  realistic  cost  associated  with  the  total  number  of  switchings  (viz.,  the  cost  of  cutting 
and  assembling  many  different  lengths  of  dielectric  material)  there  will  be  a  trade-off  between 
the  two  costs;  increasing  the  fineness  of  the  grid  on  u  (thus  increasing  the  number  of  switchings) 
without  limit  can  only  reduce  the  cost  .1  to  Jj,  that  achieved  by  the  control  Uj(z). 

One  thus  enjoys  the  consolation  that  a  candidate  solution  which  satisfies  Eqs.(3-6)  to  (3-8) 
and  (  3-16)  is  at  least  optimal  when  the  number  of  switchings  is  constrained  to  the  number  ap¬ 
pearing  in  the  solution.  To  decrease  the  cost,  the  number  of  switchings  must  be  increased; 
when  there  is  some  nonzero  cost  associated  with  the  number  of  switchings,  the  point  of  dimin¬ 
ishing  returns  is  reached  very  rapidly. 

4.  Conclusion 

The  sign  of 


at  switching  times  of  a  candidate  switching  solution  tells  whether  the  candidate  is  optimal  or 
whether  there  is  a  singular  extremal  of  lower  cost.  When  the  singular  extremal  is  better,  the 
cost  of  the  solution  may  be  decreased  by  increasing  the  fineness  of  the  grid  on  admissible  values 
of  u:  repeated  solution  of  the  necessary  conditions  (Eqs.(3-6)  to  (3-8)  and  (3-16)]  will  yield  a 
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control  function  which  has  a  cost  closer  to  the  lower  bound  Jj.  In  fact,  solution  for  successively 
finer  pa.titions  will  provide  a  measure  of  Jj  and  thus  a  measure  of  the  suboptimality  of  our 
switching  solution. 

However,  thi,s  lower  cost  is  achieved  only  at  the  expense  of  an  increased  number  of  switch¬ 
ings.  This  process  automatically  yields  a  trade-off  between  complexity  of  the  control  function 
(number  of  switchings  and  dollar  cost  of  constructing  the  indicated  structure)  and  the  cost 
achieved. 


C.  FAILURE  OF  CONVERGENCE  -  THE  SINGULAR  EXTREMAL 


In  the  previous  discussion  (Sec.  VllB),  we  assumed  that  the  numerically  derived  numerical 
solution  did  not  correspond  to  a  singular  extremal.  However,  we  found  that  in  many  eases 
(where  Eq.(7-16)  holds]  there  is  a  nearby  singular  extremal  which  provides  a  lower  cost. 

One  interesting  question  to  ask  Is,  "What  happens  when  the  numerical  solution  tries  to  track 
the  singular  extremal  exactly?  Does  the  infinitely  fast  switching  bennvlor  (approximating  Uj(7.) 
and  yielding  cost  Jj)  appear  in  the  numerical  results?"  To  answer  these  (|uestlons,  we  need 
only  examine  the  implications  of  Flq.  (7-i6);  (we  assume  in  the  following  that  relation  (7-16) 
does  hold,  i.e.,  that  our  candidate  solution  is  suboptimal). 


8u 


dz 


^  (^) 
2  'au' 


>  0 


On  this  singular  extremal,  the  multiplier  function  M  =  ail/au  is  on  a  boundary  such  that 
both  U|^  and  U|^^^  maximize  the  Hamiltonian;  if  M  were  only  slightly  larger,  the  unique  choice 
from  Eq.  (3-16)  would  be  (Here  Uj^^^  >  Uj^).  If  M  were  slightly  smaller,  the  unique 

choice  is  Uj^.  We  also  have  the  relation  dM/dz  =  0  on  this  extremal,  so  M  is  precisely  tracking 
(in  z)  the  value  which  causes  ambiguity  in  u  (see  Fig.  VlI-5). 

First,  suppose  that  we  are  "integrating"  the  state  and  costatc  equations  forward  in  z,  ob¬ 
taining  the  function  M(z)  from  x(z,  w)  and  p(z,  w)  at  each  point,  then  finding  u(z)  from  M(z).  Now 
suppose  that  due  to  finite  numerical  resolution  on  M(z),  the  control  value  Uj^^^  is  selected. 


M(l)  U-4MIU5I 


Fig.  VII-5.  M(z)  versus  z  in  the  case  of  a  singular  extremum. 
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(This  value  is,  of  course,  larger  than  the  intermediate  control  which  would  be  derived  from 
Eq.  (3-18)  and  realized  by  infinitely  fast  chattering  between  Uj^  and 
Since,  by  assumption. 


?.  2 

in  the  next  distance  dz,  d  'M/dz  will  be  >  0,  and  dM/dz  >  0,  so  M  increases  above  its  thresh¬ 
old  value  corresponding  to  the  singular  extremal.  Since  M  is  now  larger,  the  selected  by 

accident  is  now  the  unique  choice  of  Eq.(3-16);  we  have  left  the  singular  extremal.  An  identical 
argument  sliows  that  if  Uj^  had  been  selected  instead,  M  would  decrease  and  we  would  leave  the 
singular  extremal  with  unique  control  clioice  Uj^. 

Thus,  the  singular  extremal  is  unstable  with  respect  to  integration  in  the  forward  direction 
whenever  Kq.(7-lf>)  holds;  the  numerical  solution  will  very  soon  "fall  off"  this  trajectory. 

Further,  an  infinitesimal  change  in  the  initial  conditions  at  the  beginning  of  tlie  singular  tra¬ 
jectory  can  make  a  large  difference  in  the  distance  Az  before  the  solution  "falls  off"  (diverges) 
and  even  in  the  direction  in  which  the  solution  diverges.  For  this  reason,  infinitesimal  changes 
in  the  initial  state  and  costatc  will  map  into  large  changes  in  the  terminal  state  and  costate 
whenever  a  portion  of  the  trajectory  is  a  singular  extremal  (or  appears  to  be  within  the  numer¬ 
ical  accuracy  of  the  macliine). 

Interestingly  enough,  the  singular  extremal  is  not  only  unstable  with  respect  to  integration 
in  the  forward  direction,  it  is  also  unstable  when  the  equations  arc  integrated  in  the  backward 
(— z)  direction.  A  simple  repetition  of  the  above  argument  with  z  replaced  by  —  z  will  verify 
this: 

d^M  _  d^M 
dz  d(-z)‘‘ 

This  means  that  an  infinitesimal  change  in  the  terminal  state  and  costate  will  map  into  a  large 
change  in  the  initial  state  and  costatc  if  a  portion  of  the  trajectory  is  (within  measurement  ac¬ 
curacy)  a  singular  extremal. 

This  answers  questions  about  whether  fast  chattering  should  show  up  in  numerical  results 
(it  should  not,  because  independent  of  the  direction  of  integration,  the  solution  very  quickly 
diverges  from  the  singular  extremal).  Moreover,  it  points  up  an  important  limitation  of  the 
numerical  technique  described  in  Sec.  VI. 

The  basic  technique  used  there  is  an  iteration  on  the  boundary  conditions  until  all  are  solved; 
the  convergence  properties  of  this  iteration  are  based  on  a  smooth  mapping  between  initial  state 
and  costate  and  the  terminal  state  and  costate;  a  small  change  in  these  functions  at  one  end  of 
thf  line  must  map  into  a  small  change  at  the  other,  and  when  the  changes  are  sufficiently  small, 
this  mapping  is  approximately  linear.  Asa  matter  of  fact,  this  linear  mapping  shows  up  as  a 
sensitivity  matrix  in  the  Newton-Raphson  procedure  employed. 

llcwever,  when  part  of  the  trajectory  is  a  singular  extremal,  this  smoothness  of  mapping 
(and  linearity  for  small  changes)  is  lost.  The  instability  shows  up  whether  we  are  integrating 
in  the  forward  or  reverse  direction,  so  it  makes  no  difference  which  boundary  condition  we  are 
iterating  upon.  In  such  a  case,  the  procedure  of  Sec.  VI  will,  in  general,  not  converge.  Any 
scheme  which  iterates  on  boundary  conditions  is  useless  for  solving  singular  extremal  trajectories. 
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Whether  the  numerical  procedure  does  "get  onto"  a  singular  extremal  is  critically  dependent 
on  two  factors;  the  first  is  the  choice  of  the  set  of  available  u  values  {u^,  U2, .  .  .  because 

this  set  determines  the  threshold  values  of  M  where  ambiguities  in  u  occur.  A  change  in  the 
constraint  set  U  may  be  adequate  to  remove  the  singular  extremal  problem,  and  its  instabilities 
for  a  given  source  and  load.  The  second  factor  is  the  exact  values  of  the  state  and  costate  func¬ 
tions  at  the  end  of  the  line  where  iterations  begin;  a  very  slight  change  in  the  source  impedance 
may  cause  a  sufficient  change  in  these  functions  so  that  the  procedure  no  longer  "gets  onto"  the 
singular  trajectory  for  that  choice  of  {u^,  U2. .  • .  u^^} .  In  Fig.  Vll-5,  an  infinitesimal  change  in 
either  the  function  M(z)  or  ttie  threshold  value  between  different  capacity  regions  will  eliminate 
this  singular  extremal  behavior. 

Since  we  are  generally  interested  in  solving  the  design  problem  for  a  sequence  of  succes¬ 
sively  finer  grids  on  ttie  available  controls  u,  the  first  suggestion  above  is  the  more  natural 
way  to  escape  convergence  troubles  when  a  singular  extremal  shows  up.  This  suggests  the 
following  sctieme; 

During  an  iteration,  ttie  program  checks  for  simultaneous  occurrence  of  the  three  conditions: 


(a) 

M  w  a  threstiold  value  which  gives  an  ambiguity  in  u 

(7-18) 

(b) 

dz 

(7-19) 

(c) 

“  >0  . 

dz^ 

(7-20) 

WTien  all  ttiree  of  the  above  conditions  occur,  the  procedure  is  attempting  to  iterate  along 
a  singular  e.xtremal;  all  hope  of  convergence  on  future  iterations  is  lost. 

Abort  the  run  on  the  present  problem,  define  a  new  problem  with  a  different  (generally 
more  finely  partitioned)  discrete  constraint  set  {uj,  U2, . .  and  begin  iterating  with  the 

same  state  and  costate  functions  used  in  the  last  iteration.  The  new  trajectory  will  generally 
miss  the  singular  extremal  [having  provided  a  slight  change  in  IVi(z))  and  should  eventually 
converge. 

D.  RECOMMENDATIONS 

The  first  limitation  of  the  procedure  outlined  in  Sec.  VI  is  that  the  solutions  will  appear 
which  are  merely  suboptimal.  These  are  piecewise  continuous  approximations  to  the  optimal 
solution,  which  contains  infinitely  fast  chattering. 

The  recommended  procedure  is  then  to  solve  the  necessary  conditions  of  Eqs,  (3-6)  to  (3-8) 
and  (3-16)  for  a  number  of  successively  finer  grids  on  available  U,  in  each  case  indicating  the 
region  of  z  where  Eq.  (7-16)  holds.  Over  these  regions  there  is  a  nearby  singular  extremal 
having  lower  cost,  and  we  may  expect  to  obtain  improvement  by  a  yet  finer  partitioning  on  the 
range  of  available  controls.  In  the  limit  as  this  grid  becomes  dense,  we  approach  the  pre¬ 
viously  derived  global  optimal  control  Ujtz)  and  the  lower  bound  on  cost,  Jj.  However,  this 
cost  improvement  occurs  only  at  the  expense  of  increased  control  complexity;  the  person  con¬ 
structing  the  physical  realization  for  u(z)  will  most  likely  pick  a  tradeoff  point  corresponding 
to  a  very  gross  grid  on  u(z).  This  will  be  especially  true  because  the  cost  decrease  due  to 
finer  partitioning  reflects  the  amount  of  curvature  in  the  Hamiltonian;  this  curvature  is  ex¬ 
tremely  small. 
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Tlio  second  limitation  is  that  convergence  difficulties  may  appear  because  a  singular  extremal 
subarc  ap|)cars  in  the  numerical  solution  of  the  necessary  conditions  of  Kqs.(3-6)  to  (3-8)  and 
(3-16).  During  the  iteration,  such  an  occurrence  will  be  signaled  hj  the  relations  in  ICqs.(7-18), 
(7-19)  and  (7-20)  being  simultaneously  true.  In  such  a  case,  the-  run  has  little  chance  of  conver¬ 
gence  and  should  ho  aborted.  Tlie  problem  can  then  be  run  for  a  different  and  finer  partition  on 
the  available  U,  which  will  change  the  threshold  values  of  M  and  generally  cause  the  now  tra¬ 
jectory  to  miss  the  singular  extremal. 

If,  however,  duo  to  overriding  physical  considerations,  one  were  forced  to  .solve  the  design 
problem  for  the  precise  constraint  set,  source  and  load  specified,  then  the  only  reasonable  way 
to  proceed  (in  the  singular  oxtrom  a  case)  js  to  convert  to  a  numerical  scheme  which  performs 
its  iterations  in  .ho  policy  space  of  u(z)  ratl'.er  than  iterating  on  the  boundary  conditions. 

One  such  scheme  could  involve  iterating  on  the  switch  point  locations  z  where  u  changes 
value.  However,  this  type  of  schomi'  also  suffers  from  a  number  of  difficulties,  and  it  is  not 
clear  whether  such  a  procedure  would  generally  bo  an  impro%'ement  on  the  boundary-function 
iteration  method  of  .Sec.  VI. 
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VIII.  NUMERICAL  RESULTS 


A.  RESULTS  OBTAINED  ON  THE  SDS-940  COMPUTER 

1.  Introduction 

Computer  time  was  made  available  by  Bolt  Beranek  and  Newman  Inc.,  for  Implementation 
of  the  procedure  described  in  Sec.  VI.  The  computer  used  was  the  research  computer  of  the 
Interactive  Systems  Division,  an  .SD.S-940  operating  as  a  general  purpose,  time-sharing  sys¬ 
tem.  The  progiamming  language  selected  for  this  implementation  was  CAL,  Conversational 
Algebraic  Language,  an  interpreter/incremental  compiler  that  features  on-line  interactive 
program  running,  debugging  and  editing. 

Because  of  the  interactive  ability  of  the  system,  the  size  of  step  to  be  taken  for  each  suc¬ 
cessive  iteration  in  the  boundary-function  space  (the  size  of  the  vector  GRAD)  was  not  determined 
automatically.  Instead,  the  program  would  type  out  the  results  of  each  iteration  and  permit  the 
user  to  select  a  step  size  based  on  his  estimate  of  how  well  convergence  was  proceeding.  Of 
course,  in  a  later  implementation  of  the  program  operated  under  a  batch-processing  system, 
such  interaction  was  impossible  and  the  above  decision  process  was  automated  (accomplished 
by  the  program  itself). 

2.  Examples  of  Transmission  Line  .t^oblem 

The  problems  run  on  the  SD.S-940  were  all  e.xamples  of  the  transmission  line  impedance 
matching  problem  described  in  Sec.  IV.  To  make  the  problem  realistic,  the  basic  parameters 
were  taken  from  the  standard  RG-8U  coa.xial  cable: 

L(z)  =  80.16  X  10*^  henrys/ft  ,  CMIN  =  29.5  x  10*^^  farads/ft 

The  characteristic  impedance  of  this  cable  is  52.2  ohms,  and  this  was  selected  to  be  the  source 
impedance  at  z  =  -  /  for  all  frequencies  (a  resistive  source  perfectly  matched  to  the  unmodi¬ 
fied  RG-8/U  cable). 

CMAX  was  selected  to  oe  2’^CMIN,  a  condition  which  could  be  realized  by  filling  a  section 
of  the  cable  with  dielectric  l.i.ving  twice  the  dielectric  constant  of  the  standard  filler  material. 
The  frequency  band  of  Interest  was  selected  as  50  to  70  Me,  and  the  source  was  assumed  to 
have  a  uniform  available  power  density  S(u')  over  this  frequency  interval.  (This  would  be  consid¬ 
ered  a  wide-band  impedance  matching  problem  because  it  involves  a  33-1/3  percent  bandwidth 
around  the  center  frequency  of  60  Mc.) 

The  line  length  was  selected  to  correspond  to  one  wavelength  in  the  CMIN  material  at  60  Mc, 
or  10.82  ft. 

This  procedure  was  then  solved  (using  the  procedure  of  Sec.  VI)  for  a  number  of  different 
load  impedances,  ranging  from  resistances  of  8  x  R^  =  8  x  52.2  ohms  to  2.5  x  R^.  Complex 
(reactive)  load  impedances  were  not  considered. 

Figure  VIIl-l  shows  a  typical  result  oljtained  using  this  procedure.  The  source  was  con¬ 
nected  at  z  =  -f  (the  left  side  of  the  graph)  and  the  load  of  a  resistor  equal  to  2.5  R^^  connected 
at  z  0.  The  total  available  power  was  52.2  watts,  of  which  2.1  watts  was  reflected  due  to 
impedance  mismatch.  The  VSWR  (voltage  standing-wave  ratio)  between  the  source  and  load 
in  the  absence  of  an  impedance  matching  device  would  be  2.5  to  1;  this  one-wavelength  coupler 
has  reduced  It  to  1.5  to  1.  The  cost  from  the  added  term  which  is  quadratic  in  capacity  was 
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unrealistically  large  in  this  particular  example,  amounting  to  0.5  watt,  or  25  percent  of  the  cost 
due  to  reflected  power.  This  solution  should  be  run  again  with  the  quadratic  cost  term  (QCOST) 
about  ton  times  smaller  to  obtain  a  solution  which  docs  a  better  job  of  minimizing  reflected 
power  per  sc. 

The  two  functions  plotted  in  Fig.  VIII-1  are  the  capacity  multiplier  function  (repre¬ 

sented  by  x's)  and  the  capacity  C*  (z)  obtained  (the  solid  black  line).  The  center  line  represents 
the  zero  for  M^(z),  but  the  zero  for  C(z)  is  located  at  the  bottom  of  the  page  (29. 5  pf/ ft  C{z) 
59pf/ft).  The  two  fun  ons  are  overlaid  to  show  the  close  correlation  in  their  variations.  As 
indicated,  there  were  iive  available  values  for  capacity  per  unit  length,  ranging  in  equal  steps 
between  C\tlN  and  CMAX. 

We  note  that  over  intervals  of  z  when  Mj,(z)  is  sufficiently  negative,  C*  (z)  takes  on  the 
value  CMIN;  when  Mj,(z)  is  sufficiently  positive,  C*  (z)  takes  on  value  CMAX.  However,  for 
intermediate  values  of  M^(z)  (»0),  (z)  takes  on  the  three  values  available  to  it  which  arc  inter¬ 

mediate  between  CMIN  and  CMAX,  as  indicated  by  Eq.  (3-l6). 

There  is  one  particularly  important  aspect  of  this  solution;  at  the  right-ho  •  w-nd  of  the  line 
where  the  load  2.5  is  connected,  M^(z)  goes  through  large  sinusoid-like  variations  having 
successive  zero-crossings  about  1/4  electrical  wavelength  apart.  (An  electrical  wavelength 
is  times  shorter  in  the  C.MAX  material  than  in  the  CMIN  material.)  These  large  excursions 
in  M^(z)  cause  rapid  switching  of  C*(z)  between  CMIN  and  CMAX  over  successive  intervals  of 
1/4  electrical  wavelength.  If  the  quadratic  cost  term  were  made  considerably  smaller  (as  it  in 
fact  should  be)  then  the  region  of  values  giving  rise  to  intermediate  C^fz)  values  would  shrink 
to  a  very  narrow  band  around  zero.  This  would  cause  the  previously  mentioned  switchings  be¬ 
tween  CMIN  and  CMAX  to  become  sharper  (occur  over  a  smaller  distance)  so  the  C(z)  wa\ --form 
at  the  right  end  of  the  line  would  approach  a  "bang-bang"  solution,  switching  between  CMIN  and 
CMAX  over  successive  intervals  about  1/4  electrical  wavelength  apart.  This  same  behavior 
was  noted  not  only  in  the  2.5:1  match  above,  but  also  in  the  other  problems  ranging  up  to  an  S;1 
match,  and  even  in  the  waveguide  problems  described  in  the  next  section. 

To  understand  what  is  happening  in  terms  of  impedance  matching,  we  look  at  an  8:1  imped¬ 
ance  matching  method  on  the  Smith  Chart  of  Fig.  VIlI-2.  For  clarity,  we  consider  matching  at 
only  a  single  frequency,  rather  than  over  an  interval  of  frequencies.  The  source  will  be  a  re¬ 
sistor  R^,  represented  by  the  center  (A)  of  the  Smith  Chart;  the  load  will  be  a  resistor  8  R^, 
represented  by  the  point  (H).  The  available  capacities  CMIN  and  CMAX  =  2  C.MIN  correspond  to 
impedances  R^(Point  A)  and  R^/v/T  (Point  D).  Moving  down  the  line  from  load  B  toward  source 
A  is  represented  on  this  chart  by  a  circular  motion  around  a  center  specified  by  the  capacity  in 
use:  center  A  for  C  =  CMIN,  center  D  for  C  =  CMAX.  To  achieve  a  perfect  impedance  match, 
we  wish  to  describe  a  trajectory  from  B  to  A  using  the  smallest  possible  electrical  length  (rota¬ 
tion  on  chart). 

This  optimal  match  is  obtained  as  follows:  starting  at  the  load  end  (Point  B)  use  l/4  elec¬ 
trical  wavelength  of  CMIN:  this  results  in  180*  rotation  about  center  A  on  the  chart.  Now  use 
1/4  electrical  wavelength  of  CMAX,  for  180'  rotation  about  Point  D.  Repeating  this  sequence 
twice  more  results  in  spiraling  inward  from  B  to  A,  achieving  perfect  match  with  3/4  electrical 
wavelength  of  CMIN  and  the  same  for  CMAX. 

It  is  easy  to  convince  oneself  that  the  use  at  any  time  of  a  capacity  intermediate  between  CMIN 
and  CMAX  (corresponding  to  rotation  about  a  center  between  A  and  D)  will  result  in  a  slower 
inward  spiraling  and  hence  will  require  a  longer  electrical  length  to  achieve  the  same  match. 
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This  is  then  the  reason  for  the  appearance  of  "bang-bang"  solutions  switching  between  CMIN 
and  CMAX  over  successive  1/4 -electrical  wavelength  intervals  near  the  high-resistance  load: 
this  resistance  is  being  "pumped"  in  i'>ward  the  origin  by  the  same  mechanism  indicated  in 
Kig.  Vlll-2. 

Thi'i  "bang-bang"  behavior  is  only  apparent  at  the  right  end  of  the  line,  however;  at  the  left 
side  is  region  where  M  (z)  «  0,  causing  C*(z)to  take  on  intermediate  values.  If  M^(z)  is  equal 
to  zero  over  an  interval,  then  making  the  quadratic  cost  QCOSTarbitrarily  small  will  not  change 

this  behavior:  C*(z)  will  still  take  on  inter¬ 
mediate  values.  Such  behavior  appears  to  be 
the  case  in  all  examples  studied,  so  it  is  ex¬ 
pected  that,  in  general,  as  QCOST  is  made  arbi¬ 
trarily  small,  the  solution  C*(z)  will  converge 
to  a  function  that  is  "bang-bang"  over  part  of 
the  line  (like  the  successive  l/4  electrical  wave¬ 
lengths  near  the  high  resistance  load)  and  takes 
on  intermediate  values  over  the  remainder  of 
the  line. 

The  performance  obtained  in  the  8:1,  6:1 
and  1:1  resistive  matching  problems  studied  is 
shown  in  Fig.  Vlll-3.  Here  we  plot  on  the  Smith 
Chart  the  normalized  impedance  (as  a  function 
of  frequency)  seen  by  the  52.2-ohm  source  look¬ 
ing  to  the  right  into  the  coupler;  the  magnitude 
of  the  reflection  coefficient  p(ui)  is  represented 
by  the  radial  distance  from  the  origin  to  the 
impedance  curve. 

The  striking  result  is  that  the  resulting 
impedance  functions  wrap  uniformly  around  the  origin,  maintaining  almost  constant  |p(w)|  over 
the  entire  frf’quency  band.  This  behavior,  a  uniformly  good  match  at  all  frequencies  rather 
than  a  perfect  match  at  some  frequencies  and  poor  match  at  others,  was  predicted  in  Sec.  1  and 
by  the  results  of  Appendix  I. 

3.  Conclusions 

The  transmission  line  problems  studied  all  show  a  bang-bang  behavior  in  C(z)  near  the  high 
resistance  load,  switching  between  CMAX  and  CMIN  over  successive  l/4  electrical  wavelengths. 
Such  behavior  is  expected  to  be  common  in  the  case  of  severe  mismatches,  as  this  behavior 
pumps  the  impedance  level  down  toward  the  origin  of  the  Smith  Chart  (indicating  perfect  match) 
as  fast  as  possible  (in  the  minimum  electrical  length).  Over  other  regions  of  the  line’,  how¬ 
ever  (in  particular,  near  the  source  which  is  matched  to  the  CMIN  characteristic  impedance 
of  the  line),  the  function  C'*'  (z)  takes  on  intermediate  values  of  capacity  and  is  expected  to  do  so 
even  when  QCOST  is  made  arbitrarily  small.  The  resulting  matching  structure  tends  to  provide 
an  equally  good  impedance  match  at  all  frequencies  of  interest. 


Fig.  VIII-3.  Performonce  of  kxd*  tronsmiuion 
line  solufiont. 
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B.  RESULTS  USING  THE  IBM-360  COMPUTER 


1.  Introduction 

Computer  time  was  made  available  by  the  M.I.T.  Lincoln  Laboratory  Group  42  for  imple¬ 
mentation  of  the  program  described  in  Sec.  VI.  The  computer  used  was  the  IBM-i60,  Model  65, 
and  the  language  selected  for  the  implementation  was  Fortran-IV-G.  The  source  program  is 
maintained  as  a  card  deck  and  is  submitted  for  runs  under  a  batch-processing  monitor. 

Since  no  interaction  is  possible  between  the  user  and  the  program  running  in  this  mode,  the 
programmed  procedure  was  made  fully  automatic.  This  automation  was  accomplished  by  requir¬ 
ing  as  input,  besides  the  cards  specifying  the  parameters  of  the  problem,  a  sot  of  control  cards 
which  specify  the  maximum  number  of  iterations  to  be  performed,  the  maximum  step  size  in  the 
Newton-Raphson  procedure,  and  flags  indicating  whether  the  initial  guess  for  the  impedance 
presented  by  the  coupler  to  the  source  should  be  read  in  from  additional  cards,  generated  by 
the  program  itself,  or  be  set  equal  to  the  impedance  obtained  during  the  final  iteration  of  the 
previous  problem  solution. 

A  number  of  different  problems  can  be  run  in  one  session  by  stacking  as  many  control  decks 
plus  parameter  specification  decks  as  one  desires  in  the  input  deck  submitted.  The  procedure 
will  terminate  the  run  of  one  problem  and  read  in  the  next  control  deck  and  parameter  deck  to 
start  the  next  problem,  under  any  of  the  four  following  conditions: 

(a)  The  procedure  has  converged  to  the  desired  boundary  function 
impedance. 

(b)  The  procedure  has  performed  the  indicated  maximum  number  of 
iterations  on  the  problem. 

(c)  Input  parameters  are  specified  improperly;  for  example,  the  wave¬ 
guide  width  is  too  narrow  to  permit  propagation  of  real  power  at  the 
frequencies  of  interest. 

(d)  The  iteration  is  proceeding  along  a  singular  extremal  trajectory  and 
thus  has  no  hope  of  convergence. 


The  parameter  specification  deck  Includes  specification  of  the  problem  type  (waveguide  or 
transmission  line),  the  waveguide  width  if  the  problem  type  is  waveguide,  the  length  of  the 
structure,  the  Inductance  '’alue  (inductance  per  unit  length  for  transmission  line,  for  wave¬ 


guide),  the  maxi 
length  in  tran 
ation,  the 
across  this  ban> 


limum  available  values  of  the  control  parameter  (capacity  per  unit 
in  the  waveguide),  the  center  frecjuency  and  bandwidth  of  oper- 
d  load  impedances  sampled  at  eleven  equally-spaced  frequencies 
available  power  at  each  of  the  above  frequencies.  In  .addition,  an 


upper  bound  on  the  cos  Jue  to  quadratic  control  term,  QCOST,  and  the  number  of  available 
control  values,  NVALS,  are  specified.  These  available  control  values  arc  then  taken  to  be 
equally  spaced  between  the  max  and  min  available  values  specified  earlier. 

This  procedure  was  run  for  a  number  of  variations  on  a  particular  problem  to  determine 
the  change  in  cost  and  in  the  control  parameter  function  when  the  parameters  NVALS,  QCOST, 
and  length  are  changed  in  the  problem  specification.  The  results  cDtained  are  described  below. 


2.  A  Practical  Waveguide  Example 

One  problem  of  particular  interest  for  practical  application  was  studied  extensively:  the 
problem  of  obtaining  a  good  Impedance  match  between  a  narrow  slot  which  forms  one  element 
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Fig.  VIII-4.  Waveguide  solution  for  NVALS  =  2. 


Fig.  VMI-5.  Solution  for  NVALS  =  5,  (XOST  <  0. 5. 


Fig.  VII 1-6.  Solution  for  NVALS  =  20,  QCOST  <0.5. 
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of  a  phased  array  receiving  antenha,  and  a  quartz-loaded  waveguide  having  the  same  cross  sec¬ 
tion  connected  to  the  slot.  The  reeeived  power  is  then  eondueted  to  a  strip-line  reeeiver  (also 
of  the  same  cross  section)  by  Mie  quartz -loaded  waveguide. 

The  receiving  element  was  selected  as  a  narrow  slot  0.25  by  0.9  inch  to  match  the  dimen¬ 
sions  of  the  striplihe  receiver  and  thus  simplify  mechanical  construction  (eliminate  the  need  for 
‘.apered  waveguides,  etc.  ).  However,  this  slot  is  so  narrow  that  there  is  a  severe  mismatch 
in  wave  impedance  between  a  wave  in  the  waveguide  section  and  the  wave  impedance  of  the  free- 
spacc  wave  impinging  on  the  slot.  In  fact,  this  mismatch  is  so  severe  that  it  typically  results 
in  greater  than  20:1  VSWR  and  more  than  7.5  db  of  reflection  loss.  The  receiving  slot  end  of 
the  waveguide  presents  an  impedance  considerably  more  reactive  than  resistive,  and  the  imped- 
anee  mateh  is  desired  over  a  fairly  wide  (10  percent)  bandwidth  from  5  to  5.5Ge.  Such  a  diffi¬ 
cult  impedance  matching  problem  does  not  lend  itself  to  standard  techniques,  so  there  was  ex- 
eellent  incentive  to  apply  the  technique  outlined  in  Sec. VI. 

The  problem  was  run  using  the  following  constraints:  =  3-78  (quartz  die’eetrie); 

^MAX  "  ^o  dielectric);  length  =  1.15  inches  =  one  wavelength  at  center  frequency 

5.25Ge  in  quartz  dielectric.  The  available  power  was  assumed  to  be  uniform  at  eleven  equally 
spaced  frequencies  ranging  from  5  to  5.5Ge,  and  normalized  to  1  watt  at  each  frequency.  (This 
normalization  has  no  effect  on  the  resulting  (z)  solution  because  all  system  equations  and 
boundary  conditions  are  linear  in  power.) 

An  upper  limit  on  the  quadratic  cost  term  [length  x  A  ~  fg)^/2,  where  ^ 

was  selected  to  be  0.5  watt,  or  about  4.5  percent  of  the  total  available  power  of  11  watts. 
The  solution  was  th.'n  obtained  for  five  equally  spaced  available  f (z)  values  ranging  from 
to  (NVALS  =  5).  The  starting  guess  for  Ihe  impedance  presented  to  the  source  by  the 

coupler  was  generated  by  the  program;  this  impedance  would  result  if  the  waveguide  were  filled 
everywhere  with  the  dielectric. 

The  solution  required  about  10  minutes  of  run  time  on  the  IHM-360  because  the  starting 
guess  was  considerably  different  from  the  final  impedance  obtained.  The  solution  was  then 
repeated  for  NVALS  =  2,  3,  10,  15  and  20,  starting  with  the  impedance  obtained  from  the  NVALS 
=  5  solution  above.  Each  of  these  solutions  required  about  5  minutes  of  run  time,  demonstrating 
a  considerable  improvement  because  the  initial  guess  was  fairly  close  to  the  final  solution  in 
these  problems. 

The  e*(z)  solutions  obtained  for  NVALS  =  2,  5  and  20  are  showm  in  Fig.  VIII-4,  -5,  and  -6, 
respectively.  In  each  ease  we  note  the  successive  quarter-electrical-wavelcngth  seetions  of 
essentially  bang-bang  nature  (between  and  beginning  from  the  high  impedance 

souree  at  z  =  -1  .  This  is  the  same  behavior  as  that  shown  in  the  transmission  line  problems, 
and  may  be  explained  by  the  same  reasoning  (Fig.  VIII-3). 

At  the  right-hand  end  of  the  waveguide  (z  =  0),  near  the  resistive  load  having  the  charaetcr- 
istie  impedance  of  quartz -loaded  waveguide,  we  see  quite  a  different  behavior.  Here  the  ap¬ 
propriate  e*  (z)  does  not  behave  in  bang-bang  fashion,  but  appears  as  a  quantized  eontinuously- 
varying  funetion  taking  on  values  of  f  intermediate  between  and  As  the  quanti¬ 

zation  (partition  on  f)  beeomes  increasingly  finer  going  from  NVALS  -  2  to  5  to  20,  we  see 
that  this  "quantized"  solution  appears  to  be  approaching  a  continuous  solution  which  is  in  fact 
f*  j(z),  the  solution  obtained  when  f(z)  is  only  constrained  to  lie  in  the  interval 
As  In  the  case  of  any  other  numeri  al  integration  scheme,  the  burden  of  deciding  when  the  parti¬ 
tion  is  sufficiently  fine,  i.e.,  when  €(z)  has  converged  to  a  funetion  sufficiently  close  to  f'*' j(z), 
rests  upon  the  user. 
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The  performance  obtained  for  NVALS  =  S  is  as  follows:  the  VSWU  has  been  reduced  from 
20:1  to  5.4:1,  reducing  reflection  losses  from  7.5  to  2.7  db,  a  significant  improvement  over  the 
performance  without  an  Impedance -matching  coupler.  The  total  cost  obtained  is  5.5814watts, 
of  which  5.152  watts  is  true  reflected  power  and  0.4294  watt  is  due  to  the  quadratic  cost  term. 
In  all  these  problems,  the  quadratic  tost  term  took  on  about  85  percent  of  the  value  we  speci¬ 
fied  as  Its  upper  limit,  because  about  3/4  of  the  line  length  demonstrated  bang-bang  or 
'max  behavior. 


NO  OF  SALABLE  DIELECTRICS 

Fig.  VIII-7,  Behavior  of  totol  cost  versos  NVALS  for  QCOST  <0.5. 

We  stated  above  that  the  user  must  decide  when  the  partition  on  e  is  fine  enough  that  f(z) 
has  converged  sufficiently  close  to  €■>  j(z).  Another  parameter  we  may  observe  to  make  deci¬ 
sions  about  convergence  is  *be  total  cost  corresponding  to  a  partition  having  NVALS 

available  values.  ■*pjyy\L,5  is  monotone  decreasing  with  NVALS  and  converges  to  Jj,  the  total 
cost  when  e(z)  is  constrained  only  to  lie  in  the  interval  Ifjyjjpj.  '^***"^ 

Jnvals  NVALS.  First  wo  note  that  Jp^yy^Lg  is  Indeed  a  decreasing  function,  but  only 

a  3  percent  improvement  is  obtained  in  total  cost  as  the  partition  is  increased  from  two  avail¬ 
able  values  to  20.  In  fact,  most  of  that  improvement  is  obtained  in  going  from  NVALS  =  2  to 
NVALS  =  3,  because  the  addition  of  the  third  (intermediate)  value  permits  a  significant  reduc¬ 
tion  In  the  quadratic  cost  term  over  the  right-hand  quarter  of  the  line.  The  total  cost  appears 
to  have  converged  to  Jj  by  the  time  NVALS  =  20,  so  J  j  «  5. 581 2  watts.  Of  this  cost  5. 1  54  watts 
is  reflected  power  and  0.4272  watt  is  due  to  the  quadratic  cost  term. 

We  note  that  even  for  NVALS  =  2,  we  could  have  attained  this  cost  of  5.5812watts  merely 
by  permitting  infinitely  fast  switching  at  a  duty  cycle  indicated  approximately  by  the  f(z)  curve 
of  Fig.  VlIl-6.  Therefore,  the  NVALS  =  2  solution  that  was  obtained  was  suboptimal  by  about 
3  percent.  However,  this  3  percent  Improvement  which  is  possible  by  going  from  NVALS  =  2 
to  an  infinitely  fine  partition  is  certainly  not  worth  the  additional  difficulty  of  constructing 
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such  a  complex  device,  so  the  practical  implementation  would  be'done  with  NVALS  =  2  (two  dif¬ 
ferent  dielectrics),  or  at  most,  NVALS  =  3. 

Since  the  quadratic  cost  term  was  quite  noticeable  compared  to  the  true  reflected  power  cost 
(about  8  percent  in  these  solutions),  it  was  decided  to  repeat  the  solutions  for  an  upper  bound  of 
0.1  watt  on  QCOST  (five  times  smaller  than  before). 

It  was  unnecessary  to  repeat  the  solution  for  NVALS  =  2,  because  the  quadratic  cost  term 
has  no  effect  on  this  solution,  and,  in  fact,  QCOST  is  explicitly  ignored  by  the  program  when 
NVALS  =  2.  The  solution  was  repeated  for  NVALS  =  5  and  20  starting  from  the  impedance  ob¬ 
tained  in  the  original  NVALS  =  5,  QCOST  <  0.5  watt  solution.  Ilecause  of  the  change  in  the  QCOST 
bound,  A  was  five  times  smaller  and  so  was  each  threshold  value  of  the  multiplier  function  M(z). 
(A  threshold  value  is  an  M(z)  value  where  f*(z)  switches  from  one  value  to  another.) 

It  happened  that  the  function  \l(z)  was  then  almost  tangent  to  the  most  negative  threshold 
value  at  z  «-0.1  inch  (a  situation  analogous  to  that  of  Fig.  VII-5),  so  the  computer  was  almost 
iterating  along  a  singular  extremal.  Such  a  situation  caused  small  changes  in  the  impedance 
guess  at  z  -  to  map  (in  a  quite  nonlinear  manner)  into  large  changes  in  the  impedance  func¬ 
tion  obtained  at  z  -  0.  b'or  this  reason,  the  convergence  properties  of  the  iteration  were  con¬ 
siderably  worse  than  before,  ret(uiring  a  very  small  maximum  step  size  to  guarantee  convergence, 
and  about  10  minutes  of  run  time  to  obtain  this  convergence. 

This  difficulty  could  easily  be  avoided  by  further  decreasing  the  upper  bound  on  QCOST  so 
that  the  critical  threshold  value  of  .\I  is  moved  away  from  its  tangency  at  z  =  -0.1  inch  |a  local 
minimum  of  M(z)l.  In  fact,  as  QCOST  is  made  very  small,  the  convergence  properties  of  the 
iteration  will  be  identical  to  those  of  the  NVALS  =  2  solution,  which  in  this  case  were  quite 
stable  and  fast  to  converge. 

The  c(z)  solution  for  NVALS  =  5,  QCOST  <0.1  is  shown  in  Fig.  VIII-8.  Note  that  the  bang- 
bang  region  of  the  line  is  unchanged  from  that  of  Fig.  VIII-5,  except  that  the  transitions  are  five 
times  sharper  (a  transition  between  and  occurs  in  l/5  as  much  distance).  In  fact, 

as  QCOST  is  reduced  to  zero,  the  transitions  in  this  bang-bang  region  occur  in  zero  distance, 
just  as  in  the  NVALS  =  2  solution  of  Fig.  VIII-4. 

The  right-hand  end  of  the  line  shows  a  significant  change  from  Fig.  VIII-5,  in  particular, 
f(z)  is  farther  away  from  its  intermediate  value  t  for  points  z  where  t(z)  ^  f  .  This  change 
is  reasonable,  since  the  quadratic  cost  term,  A(f  —  f^)  /2,  is  a  "centralizing"  force  which 


Fig.  VIII-8.  Solution  for  NVALS  «  5,  CKOST  <0. 1. 
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Fifl.  VIII-9.  Behavior  of  total  cost  versus  NVALS  for  QCOST  <0. 1. 
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Fig.  VIII-10.  Voriotion  of  reflected  power  wi.h  NVALS  ond  QCOST. 
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tends  to  keep  f(z)  near  f  .  When  QCOST,  and  thus  A,  Is  reduced,  this  centralizing  effect  is 
also  reduced  and  the  solution  becomes  nearer  to  the  bang-bang  or  f(z)  -  or  solution 

with  the  exception  of  points  where  c  -  t  (here  the  quadratic  cost  term  is  zero,  so  the  scaling 
on  this  term  has  no  effect). 

In  fact,  for  this  particular  example,  a  further  decrease  in  A  will  only  sharpen  the  switch¬ 
ings  between  f  and  until  the  solution  converges  to  that  of  F'ig.  Vlll-4  (obtained  for 

NVALS  =  2).  However,  in  a  solution  where  M(z)  «  0  over  an  interval  of  z,  the  solution  would 
not  necessarily  converge  to  such  a  bang-bang  solution  as  QCOST  -*  0. 

For  this  example,  we  may  view  Figs.  VIII-4,  -5,  -6  as  a  progression  demonstrating  that 
f(z)  converges  to  the  continuous  function  fj(z)  as  NVALS  -•  Further,  Figs.  VIII-4,  -Sand  -5 
are  a  progression  demonstrating  that  increasing  QCOST  changes  f(z)  from  a  bang-bang  function 
to  a  smoother  function  approaching  closer  to  the  uniform  function  f(z)  =  (this  solution  is  ob¬ 
tained  in  the  limit  as  A  —  ■»). 

The  variation  of  the  total  cost  for  QCOST  <-  0.1  is  shown  in  Fig.  VIII-9  for  NVALS  =  2,  5, 
and  20.  Here  the  cost  improvement  is  only  l/3  percent  as  NVALS  is  increased  from  2  to  20, 
indicating  that  the  physical  structure  should  certainly  be  constructed  for  NVALS  =  2.  This  also 
indicates  that  our  solution  obtained  here  for  NVALS  =  2  is  only  suboptimal  by  l/3  percent. 

So  far,  we  have  only  examined 

changes  in  total  cost  =  reflected  power  +  QCOST 

as  the  number  of  available  f  values,  NVALS,  and  the  upper  bound  on  the  quadratic  cost  term 
QCOST  are  varied.  The  engineer  using  the  solution  is  only  interested  in  the  variations  of  the 
reflected  power,  so  we  examine  this  in  Fig.  VIII-10.  First,  we  note  that  for  every  condition 
examined,  the  reflected  power  does  not  vary  by  more  than  0.135  percent.  The  reflected  power 
for  NVALS  =  2  is  independent  of  QCOST,  as  is  f(z),  so  the  curves  for  QCOST  <0.5  and  QCOST 
<  0.1  meet  at  this  point.  For  the  QCOST  <0.1  curve,  there  is  no  measureable  change  in  re¬ 
flected  power  between  NVALS  =  2  and  NVALS  -  20,  again  making  it  obvious  that  only  two  die¬ 
lectrics  should  be  used  in  the  construction  of  the  coupler. 

A  feature  that  seems  surprising  at  first  is  that  for  the  QCOST  <0.5  curve,  reflected  power 
increases  with  NVALS,  that  is,  with  more  available  <  values  we  do  a  worse  matching  job. 

This  behavior  is  possible  because  we  are  not  really  minimizing  reflected  power,  but  reflected 
power  plus  quadratic  cost.  This  function  is  (and  must  be'  a  decreasing  function  of  NVALS.  The 
explanation  for  the  increasing  reflected  power  is  that  the  intermediate  values  of  e  permit  a 
noticeable  decrease  in  the  quadratic  term  (which  was  fairly  large  in  this  problem),  and  taking 
on  these  intermediate  values  of  f  reduces  the  quadratic  cost  (with  increasing  NVALS)  faster 
than  it  increases  the  reflected  power  cost.  Thus,  the  system  trades  some  reflected  power  cost 
to  attain  lower  cost  from  the  quadratic  term.  When  the  quadratic  term  is  made  smaller,  as  on 
the  QCOST  <0.1  curve,  the  tendency  toward  a  trade-off  is  greatly  reduced  and  v'e  see  no  change 
in  reflected  power  with  increasing  NVALS.  As  QCOST  -»  0,  the  plot  of  reflected  power  versus 
NVALS  must  become  a  monotone  decreasing  curve  (or  ai  least  a  monotone  nonincreasing  curve), 
so  one  CO  'Id  repeat  the  solution  for  NVALS  =  5  with  successively  smaller  upper  bounds  on 
QCOST  to  sec  if  a  reflected  pewer  cost  less  than  5. 147 watts  (that  obtained  for  NVALS  =  2)  can 
be  obtained,  or  until  one  is  satisfied  that  the  reflected  power  cost  has  converged. 
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3.  Conclusions 


The  main  conclusion  possible  from  these  examples  is  that  performance  is  essentially  inde¬ 
pendent  of  NVALS,  so  just  two  dielectrics  should  be  used  in  construction.  If  better  performance 
(leas  reflected  power)  is  desired,  one  must  relax  the  problem  constraints  by  either  permitting 
longer  structure  length  or  a  greater  range  of  f  in  the  dielectric  materials. 


Fig.VIII-11.  Solution  for  length  =  1.5  inches,  NVALS  =  5,  QCOST  <0.5. 


One  such  example  was  run.  where  the  structure  length  was  increased  from  1.15  to  1.5  inches 
(again,  NVALS  =  5  and  QCOST  <  0.5).  The  resulting  f(z)  function  is  shown  in  I'ig.  VlIl-11.  We 
note  that  the  successive  quarter-electrical-wavelength  sections  of  and  appear,  just 

as  in  all  previous  solutions,  beginning  at  the  source  (z  =  -f).  The  additional  0.35  inch  has  been 
taken  up  by  the  intermediate-value  f(z)  behavior  at  the  right-hand  end  of  the  lino.  The  cost  of 
this  solution  was  5.626  watts,  of  which  5.223  was  reflected  power  and  0.4040  was  due  to  the 
quadratic  cost  term.  This  performance  was  slightly  worse  than  that  afforded  by  the  shorter 
line  because  the  quadratic  cost  term  was  being  integrated  over  a  longer  distance,  thus  result¬ 
ing  in  a  greater  cost.  In  the  limit  as  QCOST  -•  0,  the  reflected  power  for  this  longer  line  must 
be  ■$  the  reflected  power  for  the  shorter  line. 

Coupler  performance  considerably  better  (having  a  reflected  power  less  than  1  watt)  than 
that  obtained  so  far  is  desired.  A  series  of  fut  ire  runs  is  being  planned  which  will  attempt  to 
reach  this  figure  by  loosening  up  the  constraints,  increasing  the  length  to  3  inches  and  ^yiAX 
to  16  •  (boron  nitride).  These  runs  will  be  tried  for  NVALS  =  2  and  NVALS  :  5.  but  no  sig¬ 
nificant  improvement  is  expected  for  the  increase  in  NVALS. 

C.  RECOMMENDATIONS 

Several  convergence  difficulties  were  noted  during  the  computer  runs  described  above. 

First,  quadratically  fast  convergence  was  only  observed  during  the  NVALS  -  2  run.  The  other 
runs  ended  by  exceeding  the  specified  maximum  number  of  iterations,  usually  in  a  state  of  oscil¬ 
lation  about  the  solution;  in  such  a  situation,  the  step  size  would  be  limited  by  the  specified  maxi¬ 
mum  step  size,  and  the  gradient  vector  (having  magnitude  greater  than  this  maximum  step  size) 
would  be  alternating  between  two  completely  different  vectors  on  successive  iterations.  To  shut 
down  the  extent  of  this  behavior  of  "rattling  back  and  forth"  across  the  solution,  it  was  necessary 
to  reduce  the  step  size.  This  generally  terminated  in  the  same  oscillatory  behavior,  but  in  a 
much  smaller  range  about  the  solution. 

Such  behavior  is  the  result  of  a  noticeable  nonlinearity  of  the  mapping  between  the  initial 
guess  function  for  the  impedance  and  the  resulting  impedance  function  at  z  =  0,  when  the 
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impedance  guess  is  varied  by  the  indicated  maximum  amount.  Of  course,  if  the  procedure  is 
not  attempting  to  iterate  along  a  singular  extremal,  this  mapping  is  linear  for  sufficiently  small 
variations.  Making  the  maximum  variation  of  the  impedance  guess  sufficiently  small  will  then 
cause  quadratically  fast  convergence  from  functions  near  the  solution.  However,  experience 
indicates  that  this  maximum  step  size  must  be  made  very  small  indeed  to  cause  quadratically 
fast  convergence,  wasting  a  large  number  of  iterations  (and  a  great  deal  of  time)  to  go  from  a 
distant  initial  guess  to  a  function  so  close  to  the  solution  that  oscillations  would  appear. 

A  suggested  remedy  for  this  proljlcm  is  conversion  from  the  saturating  Ncwton-llaphson 
procedure  used  now  (where  the  magnitude  of  the  variation  vector  is  truncated  wheneve^r  it 
exceeds  a  specified  limit)  to  the  segmented  Newton-llaphson  procedure  described  below. 

Kach  iteration  begins  in  the  same  way  as  the  present  iteration;  nine  scans  are  made  down 
the  line  to  find  tlic  load  corresponding  to  the  initial  guess  impedance  and  tiic  variation  of  this 
load  with  respect  to  eight  independent  variations  of  the  guess  impedance.  The  sensitivity 
matrix  is  formed,  inverted,  and  postmultiplicd  by  the  error  in  load  impedance  to  find  the  grad¬ 
ient  vector.  This  gradient  vector  indicates  the  change  in  the  guess  function  (assuniing  linear 
mapping  between  guess  and  result)  which  will  cause  the  resulting  load  to  match  the  actual  load. 

Instead  of  operating  a  saturation  function  on  the  magnitude  of  this  gradient  vector  Cl,  we 
imagine  scaling  it  by  a  scalar  r  to  produce  a  variation  of  our  guess  function  indicated  Ity  rG. 

Now  we  ask  what  value  of  the  scalar  r  minimizes  in  norm  the  distance  between  the  resulting 
load  impedance  and  the  desired  load  impedance.  First  try  r  =  1,  then  if  the  load  error  is  worse, 
try  r  -  l/2,  then  depending  on  the  previous  result,  r  =  1/4  or  r  =  3/4,  etc.  In  short,  one  per¬ 
forms  a  binary  search  over  the  scalar  variable  r  to  find  the  value  which  minimizes  the  load 
impedance  error.  Ivach  test  requires  one  scan  down  the  line,  so  if  a  half  dozen  points  are  tried, 
the  overall  time  for  one  iteration  is  increased  by  60  percent.  Once  a  "good"  choice  of  r  is 
found,  the  next  iteration  begins  from  this  point. 

This  procedure  is  expecteti  to  at  least  halve  the  number  of  iterations  required  to  obtain 
convergence  to  the  desired  load  impedance,  as  the  oscillatory  behavior  on  successive  iterations 
would  be  eliminated.  Instead  of  jumping  back  and  forth  acros  i  the  solution,  the  procedure  would 
move  to  the  solution  (to  first  order)  in  one  iteration.  Such  a  procedure  would  provide  quad¬ 
ratically  fast  convergence  from  a  much  larger  range  of  initial  guess  functions. 

The  second  difficulty  noted  was  that  it  was  possible  frr  M(z)  to  be  almost  tangent  to  one 
of  its  threshold  values,  thus  causing  the  computer  to  almc'St  track  a  singular  extremal.  This 
results  in  a  very  nonlinear  mapping  between  guess  functio.'!  •’nd  .''esulting  load  impedance,  so 
the  convergence  properties  of  the  iteration  are  very  poor.  This  problem  is  easily  avoided  by 
simply  changing  the  upper  bound  on  QCOST,  and  thus  A,  or  by  a  slight  change  in  or 

'^MIN. 

The  third  and  most  serious  difficulty  is  that  for  some  impedance  guesses  far  from  the 
solution,  the  procedure  fails  to  converge  at  all.  This  is  a  property  common  to  Newton-Uaphson 
procedures,  as  they  all  have  some  region  of  convergence;  the  initial  guess  must  be  in  this 
region  or  the  procedure  will  not  converge.  The  extent  of  this  region  is  a  direct  measure  of  the 
linearity  between  guess  and  result;  the  more  nonlinear  this  mapping,  the  smaller  is  the  region 
of  convergence. 

When  this  failure  of  convergence  occurs,  the  error  in  the  load  impedance  is  large,  but  is 
at  a  local  minimum;  to  first  order,  no  change  in  the  guess  function  improves  tliis  error,  and 
the  sensitivity  matrix  used  in  the  Newton-llaphson  procedure  becomes  almost  singular '(develops 
very  small  eigenvalues). 
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One  improvement  would  be  to  truneate  all  eigenvalues  less  than  some  threshold  magnitude 
to  zero  during  the  inversion  proeess  and  return  the  pseudo-inverse  of  this  sensitivity  matrix. 

The  iteration  proeedure  would  then  eontinue  redueing  some  components  of  the  error  veetor 
even  if  it  eould  not  reduee  all  of  them.  However,  if  the  troublesome  point  is  a  loeal  minimum 
with  respeet  to  a^  possible  variations  (and  not  just  a  saddle  point)  then  all  eigenvalues  will 
beeome  arbitrarily  small  and  the  pseudo -inverse  will  be  the  zero  matrix,  whieh  will  break  down 
the  proeedure  at  that  point. 

In  this  ease,  a  smoothing  faetor  is  needed  whieh  makes  the  mapping  more  linear  between 
guess  and  result,  and  whieh  thus  makes  the  region  of  eonvergenee  largo  enough  to  inelude  our 
initial  guess  funotion.  The  upper  bound  on  QCOST  is  just  sueh  a  faetor:  if  wo  increase  this 
bound,  A  ■*  “>,  and  the  solution  f*  (z)  -•  c  everywhere.  In  an  example  whore  an  initial  guess 
fails  to  eonverge,  we  need  merely  take  A  very  large  and  guess  the  initial  impedaneo  eorrespond- 
ing  to  f(z)  .i  r^.  Onoe  eonvergenee  has  been  obtained  for  this  A,  wo  take  a  smaller  A  and 
start  from  the  result  of  the  previous  problem.  For  a  small  enough  ohango  in  A,  the  old  solution 
will  be  within  the  range  of  oonvergonoo  of  the  new  problem.  This  proeedure  ean  obviously  be 
oontinued,  and  by  the  time  we  reaoh  the  value  of  A  we  wore  originally  oonsidoring,  wo  will 
have  a  guess  funetion  within  the  range  of  eonvergenee  of  this  problem. 

In  faot,  this  is  prooisely  the  procedure  that  was  employed  above,  i.o.,  finding  the  solution 
for  NVALS  =  5,  QCOST  <0.1  by  starting  from  the  solution  for  NVALS  =  5,  QCOST  <  0.5. 

In  short,  the  Indieated  operating  proeedure  is  to  start  with  a  rather  largo  value  for  the 
bound  on  QCOST,  obtain  eonvergenee  starting  from  some  simple  initial  guess,  then  use  the 
resulting  impedaneo  funotion  for  a  starting  guess  for  a  new  problem  having  a  oonsidorably 
smaller  bound  on  QCOST.  This  proeedure  may  bo  eontinuod  until  the  user  is  satisfied  that  the 
refleeted  power  eost  will  not  decrease  any  further.  This  proeedure  should  be  run  for  NVALS  ■  3, 
sinee  the  change  in  eost  with  inereasing  NVALS  is  extremely  small,  and  the  time  per  iteration 
is  proportional  to  the  number  of  switchings,  which  is  proportional  to  NVALS.  After  final  con¬ 
vergence  is  obtained,  one  may  start  from  this  solution  to  find  the  cost  for  NVALS  =  10,  20,  etc. 
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DC.  CONCLUSIONS  AND  SUGGESTIONS  FOR  FURTHER  RESEARCH 


A.  PRACTICAL  APPLICATIONS 

The  solutions  to  the  waveguide  problem  examined  in  Sec.  \  III  are  a  great  improvement  upon 
the  performance  of  the  receiving  slot  antenna  in  the  absence  of  an  impedance-matching  coupler. 
Further,  this  particular  problem  invoUes  such  a  severe  mismatch  over  a  wide  bandwidth  (10  per¬ 
cent)  that  standard  rule-of-thumb  impedance-matching  techniques  are  of  little  help.  The  technique 
described  in  Sec.  V  I  thus  provides  a  method  of  proceeding  to  a  practical  solution  where  other 
methods  fail.  However,  the  solutions  obtained  for  the  particular  constraints  considered  (a  length 
of  one  wavelength  and  leave  something  to  he  desired  in  terms  of  power  transmitted 

to  the  stripline  receiver  load  (almost  50  percent  of  the  available  power  is  reflected  through 
mismatch). 

For  this  reason,  future  runs  will  be  made  relaxing  the  constraints  on  length  and  U 

the  reflected  power  loss  can  be  reduced  to  10  percent  with  reasonable  values  for  length  and 
*^MAX'  expected  that  the  indicated  structure  will  be  constructed  and  tested  experimentally 

for  performance  (by  Lincoln  I.aboratory  Group  44)  in  anticipation  of  actually  using  the  design  in 
a  working  phased  array  radar. 

B.  POSSIBLE  IMPROVEMENTS 

Convergence  of  the  numerical  proceilure  should  be  improved  by  converting  from  the  saturat¬ 
ing  Newton-Raphson  iteration  technique  used  at  present  to  the  segmented  Newton-Raphson  tech¬ 
nique,  which  was  described  in  Sec.  VTII,  using  a  binary  search  along  the  gradient  vector  to  find 
the  minimum  error  point.  This  should  provide  a  considerable  improvement  in  speed  of 
convergence. 

In  addition,  the  numerical  procedure  should  recognize  when  it  is  almost  iterating  along  a 
singular  extremal  |M(z)  almost  tangent  to  one  of  its  threshold  values)  and  reduce  or  increase 
*iMAX  to  avoid  this  problem  (giving  notice  of  the  change  to  the  user,  of  course). 

Third,  the  inversion  of  the  sensitivity  matrix  in  the  Newton-Raphson  procedure  should  trun¬ 
cate  to  zero  all  eigenvalues  having  magnitudes  less  than  1  percent  of  the  maximum  eigenvalue 
and  should  return  the  pseudo-inverse  of  the  matrix.  However,  if  even  the  largest  eigenvalue  is 
less  than  some  small  number,  say  0.01,  then  the  procedure  is  caught  in  a  local  minimum  of  the 
error  function  and  the  quadratic  cost  term  (A)  should  be  increased  by  10  to  smooth  out  the  mapping 
between  guess  and  result.  This  increases  the  extent  of  the  region  of  convergence;  once  conver¬ 
gence  is  obtained.  A  should  be  reduced  to  its  original  value  and  the  iterations  continued  from 
this  point. 

In  addition,  once  a  solution  has  been  obtained,  the  following  functions  should  be  plotted  on 
the  automatic  plotter: 

1.  The  multiplier  function  M(z). 

2.  The  dielectric  solution  «(z). 

3.  The  reflection  coefficient  plw)  presented  by  the  coupler  to  the  source. 

4.  The  reflection  coefficient  pj(ai)  presented  by  the  coupler  to  the  load. 

5.  The  reflected  power  loss  versus  angle  of  incoming  plane  wave,  as 
ZSOURCE  is  varied  over  the  functions  corresponding  to  a  number 
of  different  angles  of  incidence.  (ZSOURCE  is  a  tabulated  function 
of  0)  and  0,  the  angle  of  incidence.) 
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C.  RELATED  TOPICS  FOR  INVESTIGATION 


Two  very  interesting  related  topics  present  themselves  for  investigation  at  this  point;  the 
first  is  the  question  of  the  behavior  of  the  solution  C(z)  and  corresponding  cost  J  in  the  transmis¬ 
sion  line  as  the  maximum  available  capacity,  CMAX,  is  increased  without  limit. 

Suppose  now  that  NVALS  =  2  so  C*  =  CMAX  when  M(z)  >  0.  Note  that  d^M/dz^  is  linear  in 
C  from  Eq.  (4-3 lb).  If  we  now  assume  that  there  is  some  point  z^  such  that 

dM(z  ) 

ana  from  Eq.  (4-3 lb) 

.2„ 

- =-  «  —  oC  for  large  C 

dz 


then 


M(z  +  a)  =  0 
o 


dMfz^  +  ff) 
dz 


dM(z^) 

dz 


,  dM(z  ) 

for  Ca  «  —  — - 

a  dz 


(9-2) 


when  C  is  very  large.  This  behavior 


lim  Ca  =  constant 
C— 


indicates  an  impulse  of  capacity  per  unit  length.  Such  an  impulse  of  capacity  per  unit  length 
(plotted  as  a  function  of  distance  z)  has  the  units  of  capacity  and  may  be  interpreted  as  a  fixed 
lumped  capacity  tapped  across  the  transmission  line  at  z  =  z^.  Of  course,  such  a  solution  is 
quite  attractive  physically  because  it  is  so  easy  to  construct 

The  above  rationalization  for  impulsive  C(z)  solutions  is,  of  course,  not  at  all  rigorous. 
To  develop  a  more  reasonable  argument,  we  examine  the  limits  of  the  transition  matrices  of 
Eq.  (4-11)  and  (4-23)  as  C  “  for  zf  (z^,  z^  +  a)  in  such  a  way  that  Ca  «  constant  K. 

The  result  is: 


|V|2(z^  +  c) 

- 1 

< 

N 

_ 1 

|l|^(z^  +  a) 

|ll^(z^j)  +  «^K^|v|^(Zjj)  +  n/TwK  ['/2Sj(Zjj)) 

•/SSjiZo  *  <r) 

V8Sj(z^j)  +  \fZwK\\'\^{z^) 

"/SSj^iZo  +  a) 

V8Sr(z^) 

Thus  the  magnitude  of  the  voltage  and  the  real  power  flow  are  continuous  across  the  impulse 
of  K  farads,  while  the  magnitude  squared  of  the  current  and  the  imaginary  power  undergo  finite 
step  changes  linear  in  |V|^  and  "/SSj  at  z  =  z^. 

A  similar  result  holds  for  the  costate. 
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oAz  *■  fx) 

’  1  o 

pJz  *■  rr) 

.  ' 


P, (••«„)  ♦  «‘^K^P2 

'’z'-'o' 

PjiZjj)  -  v'^wu’Kp^ 

P4<V  J 


(9-4) 


Thus,  p^  P4  continuous  while  pj  and  p^  undergo  step  changes  linear'  in  and  p^.  .Since 
we  are  considering  the  case  of  WAU.S  -  2.  there  is  no  reason  for  a  quadratic  cost  term;  this 
means  that,  as  shown  in  Sec.  Ill,  the  cost  .1  (which  is  now  exclusively  reflected  power)  may  he 
evaluated  as  p'(z)  x(z).  Evaluating  the  cost  on  both  sides  of  the  impulse  using  the  above  state 
equation  and  analogous  costate  equations  (9-3)  and  (9-4), 


.l(z  +  (t)  -  .I(z)  -  p'(z  (t)  x(z  (t)  —  p’(z  )  x(r;  )  0 

o  '  o  o  '00 


(9-S) 


which  was  ccrrtainly  to  be  expected  if  the  solution  is  to  be  physir-ally  meaningful.  A  further 
evaluation  establishes  that  the  Hamiltonian  is  also  continuous  across  the  impulse  of  capacity, 
so  the  impulsive  solution  seems  to  be  a  perfectly  proper  and  I'casonable  solution  to  the  problem. 
To  really  develop  this  idea,  two  phases  of  work  would  be  neces.sary;  first,  one  must  foi-- 
mulate  a  distributed  parameter  maximum  principle  for  impulsive  controls,  then  apply  this  prin¬ 
ciple  to  the  transmission  line  problem  being  considered  as  CMA.X  -  to  demonstrate  that  im- 
[lulses  of  capacity  per  unit  length,  i.e.,  lumped  capacitors,  are  a  valid  solution  to  the  problem. 
Then  one  must  modify  the  procedure  given  in  .Sec.  \  1.  in  particular  causing  the  scan  down  the 
line  to  find  the  points  z.  where  .\l(Zj)  ’  0.  d.\l(z.)/dz  0.  and  there  ii  serting  lumped  capacitors 
of  value 


-2(dM  'dz)  (z.) 
(V.X'  (d^.\l/dz^) 


•Lrfl  '  ’‘3'’2’ 


(9-6) 


with  the  x's  and  p's  evaluated  at  z.  —  <r.  Then  the  appropriate  stc))  change  must  be  made  in  the 
state  and  costate  (I'lqs,  (9-3)  and  (9-4)  and  the  scan  down  the  line  continued  to  find  tbe  next  point 
where 

d.M(z.) 

M(z.)  0  ,  — r—^  '■  9 

1  dz 

Such  a  c  hange  is  not  a  very  di-astic  modification  to  the  program,  and  seems  well  worthwhile 
in  terms  of  transmis'sion  line  couplers  which  are  extremely  easy  to  c-onstruct. 

The  other  topic  which  seems  well  worth  in .  estigating  at  this  |)oint  is  the  design  of  lossless 
ladder  impedance  matching  networks  using  a  "discrete-time"  distributed  narameter  maximum 
principle.  Consider  the  lumped  1,-C  ladder  network  shown  in  Kig.  I\-1.  with  complex  voltages 
and  currents  defined  over  the  rung  index  as  V(n)  and  l(n).  The  ladder  is  assumed  to  have  a  fixed 
number  of  rungs.  \. 
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X  (w)  |I  -  A  ,(u),  L  ,))  X  A^-) 

n  '  n-r  '  n-1  n-1  '  n-1 


Wo  must  now  formulate  a  clist  rote-timc  distributed  parameter  maximum  principle  and 
verify  the  equivalent  of  directional  convexity  for  the  system  Mq.  (9-9).  Assuming  that  this  prin¬ 
ciple  parallels  the  ordinary  difference  equation  maximum  pi’inciple,  there  will  exist  (as  a  nec¬ 
essary  condition)  a  co-state  p  (w)  which  forms  a  scalar  Hamiltonian  as  follows: 

n 


ll(x  ,  p  ,  C  ,  n) 

n  '  n+r  n  n 


|)'  .  .(ai)  A  (w,  I.  . 
'  n^r  n  n 


C  )  X  ( w)  da^ 
n  n 


(9-10) 


We  note  that  this  Hamiltonian  is  quadratic  both  in  1.^  and  C'^^,  so  the  requirement  of  minimiz¬ 
ing  the  Hamiltonian  with  respect  to  I-  and  ('  will  alwi  v»  yield  unique  solutions  for  L  and  Cl  . 

2  n  n 

(The  coefficients  of  1.^  and  in  the  Hamiltonian  may  be  interpreted  as  magnitudes  squared  of 

a  power  flow,  and  hence  are  real  and  positive.)  In  this  case  the  addition  of  a  term  quadratic  in 

the  control  parameters  to  avoid  the  complexities  of  singular  extrema  is  unnecessary. 

It  appears  that  both  the  necessary  theory  and  the  details  of  an  iterative  numei-ical  scheme 
for  the  necessary  conditions  (perhaps  even  sufficient  conditions  would  he  available  for  this 
problem)  are  a  feasible  future  development. 

If  such  a  solution  is  feasible,  it  would  be  of  considerable  practical  utility  because  such  a 
matching  netu  ork  would  be  quite  easy  to  construct.  The  next  obvious  extension  would  bg  a 
conversion  to  filter  networks  rather  than  impedance  matching  networks. 
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APPENDIX  A 

OPTIMAL  CHOICE  OF  REFLECTION  COEFFICIENT 


Civcn  the  constraint 


^  In  t/|p(a')|  do)  C  M 

*  O 

we  wish  to  demonstrate  that  a  local  minimum  of  the  functional 


■'=,1 


|p(u,’)l  do? 


is  obtained  by  the  strategy 

=  I 


oxpl- M/W)  ,  uicW 
1  . 


Consider  a  strong  first-order  variation  about  this  strategy  of  th('  following  form: 
f  expl-  K/a)  .  ccta 


|p(w)|  = 


exp  (- 


(a  +  1)1  M  -  K\V 


I  . 


expl- M/U'l  .  u.’€{\V  — a  — b) 

1  .  w/W 


where 


acW  ,  beW  ,  acb  =  0 


It  is  simple  to  confirm  that  this  variation  satisfies  the  integral  constraint,  i.e., 


^  lnl/|p(w 


)|  du.-  =  M 


and  the  perturbed  cost  is 


\  |p(u!)|^  dw  ^  a  expl-  2K  a)  +  b  exp(-|(a  f  b)  M  -  2KWl/\V} 

J\V 

+  (W  -  a  -  b)  expl-  2M/\V1 

Differentiating  with  respect  to  K  and  setting  the  derivative  to  zero  will  give  a  stationary 
value  for  the  perturbed  cost.  This  results  in 

K  -  aM/\V 

which  means  that  the  only  stationary  pi-rturhation  in  ipio;)!  is  the  zero  perturbation.  To  show 
that  the  unperturbed  'p(w)'  characteristic  gives  not  only  a  stationary  cost  but  a  minimal  cost, 
we  differentiate  again  with  respect  to  K 


3  .1  d  ,  2M,  .  . 

3K  K  =  (aM/\\) 


Therefore  the  unperturbe  l  sti-ategy  'p(w)'  -  e.\p(-  .M/W),  jjcW  is  indeed  a  strong  local  minimum. 


» 

r. 


k 


r 

IV' 
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APPENDIX  B 

FLOW  CHARTS  FOR  ITERATIVE  NUMERICAL  PROCEDURE 
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MAIN  HOCKAM  -  lUHATION  IN  FUNCTION  SPACl 


Fig.  B-1.  Iferation  in  Tchebycheff  space. 
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MAIN  MOGKAiA  -  SCAN  DOWN  LINE 


|l-<»-11lll| 


SIAKT 


Fig.  B-2.  Main  icon  downline. 
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«IIU«N 


Fifl.  B-3.  Subfooh'ne  DISH. 
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SueiiOui:N[  swriCH  CMUil.  lOC,  Dist,  rjtwvAi.. 


"[5-«Miir^ 


ICIU«N 


Fig.  B-4.  Subroutine  SWITCH. 
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|5-4;-IiU0l 

SUtKOUIIN(  SUf  (DISI2,  CMUII) 


1 


USE  TRANSITION  MATRICES  (EO.  4.11  AND  4.23) 

TO  RRORAGATE  STATE  AND  COSTATE  (XOLD,  POLD) 
DOWN  LINE  (TO  XNEW,  PNEW)  bV  DISTANCE  DIST  2 

RECOMPUTE  multiplier 
(CMULT,  EO.  4.28)  AND 

Function 

ITS  SLOfE  (EO.  4.3IA) 

7 


UIUKN 


SUiKOUIlNI  VALUC  (CMUII,  NtWVAl) 


Fig.  B*5.  Subroutine  STEP  ond  VALUE. 


I  I 


SUWtOUIlNt  MATINV  (M,  «OW5,  CCXS) 


OriOINAl  MAIDIX  Mll.JI  fO«  I  •  I  10  «OW5 
fO«  J  >  I  lO  kOWV 

INVtm  »Y  SUCCtSSIVt  GAUSS  «ft>UCTIONS, 
KtOUCING  lY  MAX  AtSOlUU  tltMCNT  AVAILAILt. 

INVWSI  MAIIIX  IS  Iffl  IN  M(I,J),  I  -  I 

lo  «ows,  J  HOWS  •  I  lo  cois  ■  2  •  rows 


Fig.  B-6.  Subroutine  MATINV. 
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APPENDIX  C 

LISTING  OF  FORTRAN-IV  SOURCE  PROGRAM 
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w  o  o  o 


C  HAIN  PROGRAM 

LOGICAL  MAVGOEf  FORFLG 
LOGICAL  SELFST.CONIINt  TYPOUT.PUNCH 
INTEGER  OLOVAL 
INTEGER  PASS 

REAL  LtLENGTH.LOCtKAXSTPtIMAGPT  .HAXJHP.MAG 
CCMPLEX  ZSOURC(ll).ZLOAOtll>*TSTRTIlll.ZESTIlll.JJ 
DIMENSION  QUAOIA.ill.TCHEBIli.A)  .  SNSMAT 1 8 . 1 6 ) . TLOAO I  6  )  ,  TUUFF ( d ) • 
1GRA0I8),PWR0EN( 11 ltSECVAR(200) 

COMMON  NVALSt THRESHtLtC (PI t SORT? tFORFLGt SLOPE .CURV.MAXSTP, 

IFREOI 111 .KNEW ( 11,4I.PNEW(11.4).X0LD(I1.AI,P0L0( 1 1 . A  I . RCHAR ( 1 1 ) , 
2NSTEPSt01ST2,CMIN,CMAX»MAVG0E<WlDTH.CAPI200l tSWCHPT(200) 

NAMEL I  ST/POL Y/CUAO«TCHEB 

NAMF.LIST/CNTRL/SELFST(C0NT1N«TYP0UT.MAXJMP.  ithax.punch 
NAMELIST/PARAMS/  WAVGOE « NVALS t LENGTH.MAXSTP tL t CMI N . CM AX, UCOS T, 
1FREQM,BNOHTH,PNROEN,2SOURC,ZLOAD,VAR  /N I DTH 1/ W I OTH 
NAMELIST  /ZSTART/ZSTRT 
CALL  MSGOP 


SETUP  CONSTANTS 
PI-3.1A1593 
SCRT2-SQRT(2.0) 

JJ-(0.0  •  1.0) 

NERKS«0 
C 
C 

76  REAr(S,P0LY,ERR«90) 
tiRlTEI6,123) 

123  FCRMATIM  QUADRATURE  MATRIX  -  *,/) 

DC  79  I«l,4 

79  WRITE(6.130)IQUAD( I.J) ,J«1.11) 

130  FCRMATIlX.llFll.A) 

NRITE(6.103) 

103  FCKMATI/,*  SAMPLED  TCHEBYCHEFF  FUNCTIONS  ■  •,/) 

DO  80  I>1,11 

60  MRITE(6.131HTCHEBn,J),J«l,A) 

131  FORMAT! lX.AFll.A) 

C 

C  BIGGEST  LOOP  OF  ALL 

70  hRiTEIB.lAA) 

166  FQRMATI'l  ••••••♦•  BEGIN  NEW  PROBLEM  ••••••••  •,//) 

REAOI5.CNTRL,ERR>9l,ENO-a6) 

C  READ  IN  PARAMCTERS 

77  HEA0(S,PARAMS,ERR«92) 

11  hRITEI6,CNTRL) 

WRITEI6,PARAMSI 
IF(.N0T.NAVG0EIG0  TO  63 

78  REAOIS.WIOTHI) 

NMITE(6.WI0THn 

IF!L*CMIN*!2.*WI0TH*1FRECM-BNDWTH/2.))**2.GT.1.)G0  TU  03 
$5  hriITEI6.139l 

139  FORMAT!///, *  WAVEGUIDE  IS  TOO  NARROW  TO  PROPAGATE  REAL  POWER  AT 
lYCUR  LOWEST  FREQUENCY.  CHOUSE  WIDER  WAVEGUIDE  TO  GET  AHOVt 
2CLT0FF. »,//•  PROBLLM  IS  ABORTED.  WILL  TRY  NEXT  PROBLEM....*,/// 
GC  TO  70 
C 
C 

83  DC  36  1-1,11 

56  FREQ! I l-FRFQM»BNOWTH* ! ! I-6.0I /lO. ) 

C 

C 
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ooo  r>  o  o  r>  o  ^  onoor>  o  r)  o  r> 


IFISELFST)G0  TO  72 
IFICONTIN)GO  TO  73 
66  RE*OI5,ZST*<tT,ERR«93) 

73  CCNTINGE 

hRITEI6t2ST*RT) 

0ELT*«8.*QC0ST/ILENGTH*ICM*X-CHINt**2l 
TFRESH»DELTA*ICMAX-CMIN)*INVRLS-2. »/l2.*NVALS-2.) 
ThRSH2«.D5*OCOST/<CMAX-CMIN) 
hRITE(6t 1A3)0ELTA,THRESH,THRSH2 
IA3  FORMAT*/,*  DELTA  ■  *,£15. 5,*  THRESH  «  *,E15.5, 

I*  TMRSH2  -  *,£15.5,/) 

ITER«0 
DC  30  I>UA 
TLOAOt  n«o. 

TLOAOl I*A)«0. 

DC  30 

TLOAOt I )-TLOAOI  I  )  >  QUAD  ( I  ,  J ) •RE AL i 2L0A0IJ ) ) 

30  TLOADI  I«A)«TL0A0n>4)40UA0(  I  •  J  )«AIMACI /LOAD  I  J  )  ) 
bRITE(6,U0)(TL0A0(  I  ),  I«1,RI 
110  FOKMATI*  TLOAO  '  '.SEIS.A) 


VERY  BIG  LOOP 
51  OC  31  1-1,11 
31  ZESTII )-2STRT(  n 
MRITEIA, 113) ITER 

113  FCRMATCl  ITERATION  -  *,1^) 

CALL  TIMHMSI IH,IM, IS) 

IFITYP0UT)HRITEI6,122)IH,IM,IS 

122  FORMAT!*  HOURS  -  * ,  I A , 3X , *M I NUTES  •  *  ,  I A ,  3X, *  SECONDS  *  *,IA) 


PASS  ZERO 


PASS-0 

IF(TYP0UT)HRITEI6,109) 

09  FORMAT!///,*  PASS  ZERO*,///) 

GC  TO  60 

32  OC  33  I-1,A 
TBUFF! I )»0. 

TEUFF! I*A)»0. 

OC  33  J«l,ll 

TBUFF! I ) -TBUFF!  I  )tOUAO! I , J ) •XNEW! J,A ) /  !  SORT 2*XNEH! J , 2 ) ) 

33  TBUFF! ItA)-TBUFF! ItA)«CUAO!T, J)*XNEW!  J,3)/l S0RT2*XNEW( J,2) ) 
NKITE!6,111)!T0UFF!I ),I-1,8) 

111  FORMAT!///,*  TBUFF  •  *,8F11.A,//) 


PASS  ONE 


DC  3A  1-1,11 

3A  ZEST!I )«ZSTRT!  I  )»VAR*TCHEB! 1,1) 
IF  !TYP0UT)HRITE!6, 112) 

112  FCRMAT!///,*  PASS  ONE*,///) 

GC  TO  60 


o  o  o  o  o 


35  DC  36  l-l.ll 

36  ZESTI1I-ZSTRT( n«VAR6TCHEB(I.2l 
lF(tVP0UnMRITE(6»116) 

116  FCRHAT(///«'  PASS  TUO*t///) 

GC  TO  60 


PASS  THREE 


37  DC  38  1-1.11 

38  ZEST(1I-ZSTRT( 1 l♦VAR•TCHE8( I .31 
IF(TYP0UT)MRITEI6.115) 

115  FCRMATl///.*  PASS  THREE*.///) 

GC  TO  60 

C 

C 

C  PASS  FOUR 
C 

39  OC  60  1-1.11 

60  ZESTin-ZSTRTl  1)^VAR*TCHEBI1.6) 
1F(TYP0UTIHR1TEI6.116) 

116  FORMAT!///. •  PASS  FOUR*.///) 

GC  TO  60 

61  OC  62  1-1.11 

62  ZESTd )-ZSTRT( 1 l♦JJAVARATCHEB(  1 . 1 ) 
C 

C  PASS  FIVE 

1FITYP0UTIHR1TE(6.117I 

117  FCRMATI///.'  PASS  FIVE*.///) 

GO  TO  60 

63  00  66  1-1.11 

66  ZESTI 1 )-ZSTRT( 1 )^JJ6VAR6TCHEB( 1.2) 
C 

C  PASS  SIX 

1F(TVP0UTIMR1TEI6.118) 

118  FCRMATI///.*  PASS  SIX*.///) 

GO  TO  60 

65  OC  66  1-1.11 

66  ZESTI  1  l-ZSTRTI  1 )♦ JJ6VAR6TCHEBI 1 .3) 
C 

C  PASS  SEVEN 

1FITVP0UTIMR1TEI6. 119) 

119  PORMATI///.*  PASS  SEVEN*.///) 

GO  TO  60 

67  00  68  1-1. 11 

68  ZESTI  1  l-ZSTRTI  1 I♦JJWVAR•TCHEBI I .6) 
C 

C  PASS  EIGHT 

1FITYP0UTIMR1TEI6.120) 

120  FCRMATI///.*  PASS  EIGHT*.///) 

GO  TO  60 

C 

C 

C  MAIN  SCAN  OOMN  THE  LINE 
C 

60  LOC— LENGTH 

S)iCHPTIl)-LOC 

NSTEPS-1 

lCNT-0 

FCRFLG-^TRUE. 

FCRSGN-1.0 

CCSTd-0. 
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CMULT>0. 

0IST2>0. 

DC  1  I>1.11 

PCLO( 1,1 )>!./( 8. *RE Alt :SOURC( 111) 

PCLO(  I  •2I>(REALI 2S0URC( I  I  I **2» A I  MAG C 2S0URC I  I  I l•*2l/ 

1(8.*REALI /SOURCI  I  I  I  I 

PCLOl  I  ,  3I>S0RT2*A1MAGI  2SOURC(  1  I  I  / ( 8. *RC Ai.  t  2  SOURC I  I  I  I  I 
PCL0(I,4)«-S0RT2/e. 

XQLOl I ,4I>SQRT2*8.*PWRDENI I )*REAL(  2SQURC( I  I  I  *RF AL I  2 EST I  I  I  I / 

1( (REAL(2S0URC( I l♦^EST( 1  I  I  I **2« I Al MAGI  2 SOUHC  I  I  l«2ESTIIIII**2) 

XCLOl I • 3I-XOLO( I •A)*AIMAG(2ESTn 1 1 /REAL  I 2EST I  I  I  I 

XCLDI I •2l«XOLO( 1,41/1 S0RT2*REAL  I  2EST  I  I  I  I  I 

XCLOl I , l)«XOLO 1 1 ,2)*IREALI 2EST(  1  I  I  •*2» AI MAG  I 2E ST  I  111**21 

1  CMULT>CMULT*2.0*S0KT2*PI*FRE0I  I  I  * (  XOLO I  I , 3  I *POLDI 1,21* 

IXOLOl I , 1 l*POLU( 1 ,311 

DC  24  1>1,11 
DO  24  J>1,4 
PNEW(l,JI>POLO(I,JI 
24  XAEW(l,JI>XOLO(l,JI 

CALL  VALUEICMULT,NEhVALI 

C*CM|N*( (CMAX-CMINI /INVALS-ll )*NEWVAL 

CAP! 1 l-C 

OLOVALaNEWVAL 

DC  20  I >1,11 

IFIMAVGOEI  GO  TO  19 

T1«C 

CC  TO  20 

19  T1«C-1./IL*(2.*WIOTH*FREOII ))**2I 

20  RCHAR(1I*SQRT(L/T1) 

PI.RREF»0. 

HAM1LT>0. 

PkRNET*0. 

PliRAVL-<0. 

DO  3  I  >1,11 
OC  2  J>1,4 

2  PliRREF>PWRREF*POLO(  I ,  J)*XOLO(  I  ,J) 

Pl,RAVL>PWRAVL*PWROEN(  I  I 
PhRNET>PWRNET*XOLO( I  ,4 1 / 1 2. •S0RT2 1 

3  HAMILT>HAM1LT* IC*( XOLOl I , 3 1 *POLO « I  , 2 1 ♦ XOL 01  1 , 1 )*PDLDI  1,3)1- 
1L*(X0L0( I ,3)«POLO( I  ,1 )*XOL0( I ,2)*P0L0(  1,31  I )*FRF0( I  I 

VSMR>( l.*SQRT(PMRREF/PMRAVL) )/( 1 .-SORT ( PMHREF / PWRAVL I  I 
HAMILT»2.*SQRT2*PI*HAMILT-0ELTA/2.*«C-CMAX/2.-CMIN/2. 1**2 
AITEN>10.*ALOG10(PMRAVL/PMRNET) 

IF((. NOT. TYPOUTI.ANO. (PASS. NE.OIIGO  TO  4 
WRITE  I  6, 100)PWRAVL,PWRREF,PWRNET,VSWR,HAMILT,ATTCN 
100  FORMAT!/, •  AVAIL.  PWR  >',F6.3,»  WATTS  REFLECTED  PWR.  >  •, 
1F6.3,3X,'WATTS',3X,«NET  POWER  OELIVEREO  > ' , F6. 3, 3X , • WATTS •  , // , 

2*  EFFECTIVE  VSWR  >  «,F10.4,>  TO  1 .  *  ,  5X , *  HAM  I L TONI  AN  >  ',013.4, 3X, 
3'WATTS  PER  UNIT  01  STANCE ',//, '  NET  ATTENUATION  >  »,H0.4,'  OB',/1 
WRITE(6, 1081 
OC  18  I>1,11 

REALPT>XNEW( I,4l/ISaRT2«XNEW(l  ,21) 

IHAGPT>XNEW( I , 3) /( S0RT2*XNEW(  1,2)1 
18  WRITE(6,106I  F REQ ( I  ) , RE ALPT ,  I MAGPT 
LOGICAL  DUMP 

DOUBLE  PRECISION  0VAL,0INT 

NAMELIST/BARF/LOC,()IST,OIST2,CMULT,SLOPE,CURV,FORFLG,OLOVAL, 
1NEWVAL,NSTEPS, ICNT,OVAL,OINT,ERR,TOL 
NAMELI ST/FRAB/OLOVAL,NEWVAL 

NAMELI  ST/8LAH/CMULT,SLOPE,CURV,OIST,LOC,NSTFPS,OLDVAL,NEWVAL, 
10IST2,ICNT,OVAL,OINT,ERK,TOL,FORFLG 
C 
C 
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r>  r>  r> 


MAIN  LOOP 

DECIDE  HOW  FAR  TO  STEP 
A  CALL  OISTKCMULT.LOCtOlST.LENCTH) 

IF(CAP(NSTEPSI.EQ.CAPfNSTEPS-inWRITE(6.BLAH) 
IF(ABS(0ISTI-10.**(*SI IS. 5.6 
C  SWITCH  TO  DIFFERENT  CAPACITY  INTERVAL 

5  CALL  SWITCHICMULT.LOC.DIST.NEWVALI 
IFILOC.EQ.O.IGO  TO  ID 
FCRSGN-1.0 

IFINSTEPS.LT.200IG0  TO  S9 
NRITEI6. 1471 

147  FCRMATI/t*  YOU  HAVE  MORE  THAN  200  SWITCHINGS  IN  THE  LINE  AND 
lARRAYS  ARE  OVERFLOWING. •./•  *  DECLARE  A  COARSER  GRID  ON  AVAILABLE 
2CAPACITIES. './/.*  PROBLEM  IS  ABORTED.  WILL  READ  NEXT  PROBLEM 1 

GO  TO  70 
S9  VAR2-0. 

DUMP-. FALSE. 

DC  16  I-l.ll 

16  VAR2>VAR2«(FRE0llt**3*(XNEW(I.ll*PNEWI 1 .3WXNEWI I.3I*PNEW(  1.2)11 
SECVARf  NSTEPSI— 16*S0RT2*PI»*3*L*VAR2 

C0STQ«C0STQ4(SWCHPTINSTEPS)-SWCHPT(NSTEPS-in •DELTA/2.* 
llCAP(NSTEPS-ll-CMAX/2.-CMIN/2.l**2 
OLOVAL-NEWVAL 
ICNT-0 

1F(ABS(SL0PEI-.01*THRESH/LENGTHIS7.S7.69 
69  IF( (THRESH. NE.0).0R.(ABS(SL0PE).GT.THRSH2/LENGTHIIG0  TO  6 

57  WR1TE(6.140I 

140  FCRMAT(///«<  YOU  ARE  ATTEMPTING  TO  ITERATE  ALONG  A  SINGULAR 

lEXTREMAL  TRA.JECTORY.it  IS  UNSTABLE  AND  YOU  HAVE  NO  HOPE  OF 

2C0NVERGENCE.  CHANGE  YOUR  CAPACITY  CONSTRAINT  SET  AND  TRY  AGAIN.'. 
3//.'  PROBLEM  IS  ABORTED.  WILL  TRY  NEXT  PROBLEM....  './//I 
GC  TO  70 

6  LCC-LOC*OIST 
0IST2«0IST2*0IST 

C  PROPAGATE  STATE  AND  COSTATE  DOWN  LINE  BY  DISTANCE  'OIST* 

9  CALL  STEP(0IST2.CMULTI 
ICNT-ICNT*! 

IF( ICNT.GT.20I0UMP«.TRUE. 

IF(ICNT.GT.30IWRITEI6«148I 

148  FORMAT! /*•  HANGUP  •  MORE  THAN  30  ITERATIONS  ON  A  SINGLE 
ISWITCHPOINT.'./I 

IF( ICNT.GT«30)G0  TO  70 
IFITHRESH.NE.OIGO  TC  88 

IF(  (ABS(CMULT).LE.THKSH2).AN0. (ABS(CMULT*SL0PE*LENCTHI.LE.THRSH2n 
IGC  TO  5 

IF(  (-CMULT^FORSGN/SLOPEI.LT.OIGO  to  58 
I F  ( -CMUL  T*FORS GN/  SLOPE- lOE-5 15.5.58 
88  OV AL> 1. 000* (-1.000*NVALS/2. 0001 •( 1.000* (0BLE(CPULTI«10.0D-5* 

1A8S(  SLOPE) l/THRE SHI 
IF  IOVAL.lt. l.OOOIGO  TO  58 
IFI0VAL.GT.(NVALS-.500))G0  to  58 
OINT>IOINT(OVALI 
ERR>OVAL-OINT 

TCL> I NVALS-2.0l*10.0E-5*ABSf SLOPE  I /THRESH 
IFIERR.LE.TOLIGO  TO  5 

58  IFI0UMP)WRITEI6.BARF) 

CALL  VALUEICMULT.NEWVAL) 

C  CHECK  FOR  ENO  OF  LINE 

IFIILOC  .EQ.  0.1  .AND.  (NEWVAL  .EQ.  OLOVALII  GO  TO  10 
IFI IABSINEUVAL-OLOVAL)-1)  7.7.8 
C  OVERSHOT  TWO  OR  MORE  SWI TCHPOINTS-BACK  UP 

8  OIST-OIST/2. 
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LCC*LOC-DIST 
0IST2»PIST2-0IST 
GO  TO  9 

C  IF  NFVtVAL  .NE.  OLOVALi  OVERSHOT  SHI  TCHPT-ST  ART  BACK 

7  IFINEWVAL.EQ.OLOVADGO  TO  89 
FCRFLG-.NOT.FOkFLG 
FCRSGN=-FORSGN 
89  IF(CUMP)MRITE(6.FRAR) 

CLOVAL'NEHVAL 
GO  TO  A 
C 
c 

C  FINISHED  WITH  A  SCAN  DOWN  THE  LINE 
C 

10  IF( .NOT.TVPOUT IGO  TO  75 
V«H  ITE(6.  1071 

107  FORMAT!/, •  LOCATION  CAPACITY*,/) 

OC  12  l-l,NSTEPS 

12  hRITE(6,105ISWCHPTin,CAP(I> 

105  FCRMATI  lX,F10.6,eX,E15.<il 
HAMILT>0. 

00  7A  1^1,11 

7A  HAMILT«HAMILT«(C*IXNEWfl,3l«PNEWII,2)«XNEHl 1 , 1 ) APNE W I  1 , 3 li¬ 

lt*  I XNEM  (  1 ,  3 1  APKEV,  I  1 , 1 1  4  XNE  W I  I  ,2)*PNEW(  1,3) I )*FRE0l  1 ) 
HAMILT«2.*S0RT2*PI*HAMlLT-nELTA/2.*IC-CMAX/?.-CMIN/2. )**2 
WRITE(6, I2eiHAHlLT 

128  FORMAT!*  HAMILTONIAN  >  *,610.5) 
hR ITEI6,  lOei 

loe  FCRMATI3X,*FRECUENCY*,10X,*REAL  PT . * • SX, *  1  MAG.  PT.*,/) 

OC  17  1-1,11 

REALPT-XNEM!  I  ,  A I  / !  SOR  T2*XNEHn  ,2 )  ) 

IMAGPT-XNEHI l,3l/(SQRT2*XNEW! 1,2)) 

17  WKITE!6,106)  FREQ! 1 1 ,REALPT, IMAGPT 

106  FORMAT! 1X,E  12. A, eX,F 10.5, 6X,F 10.5) 

75  COSTQ-COSTQ-SWCHPT INSTEPS) •0ELTA/2.*ICAP|NSTEPS)-CMAX/2.-CMIN/2. ) 

1**2 

IF!TYPQUT.OR.!PASS.EQ.O) )MRITEI6,12A)CCSTQ 
12A  FORMAT!/,*  COST  ASSOCIATED  WITH  QUADRATIC  CAPACITY  TERM  ■  *, 

IFF. A,*  WATTS*,/) 

IF(PASS)63,62,63 

63  OC  65  1-1,4 
SNSMATII ,PASS)-0. 

SNSMATI l44,PASS)-0. 

DC  61  J-1,11 

SNSMATI I,PASS)>SNSMAT<  I,PASS)*QUAOII,J)*XNEMIJ,A)/ISORT2*XNEMI  J,2) 
1  ) 

61  SNSMATII44,PASS)«SNSMAT<I44,PASS)*QUA0II,J)*XNEWI  J,3)/IS0RT2* 
IXKEHf J,2) I 

65  CONTINUE 

IF! TYPQUTIMRITE<6,121) <SNSMATII,PASS) ,I-I,B) 

121  FCRMATI//*  test  -  *,8F11.5,//) 

OC  64  1-1,4 

SNSMATI I,PASS)-ISNSMAT<I,PASS)-T6UFFI I ) )/VAK 

64  SNSMATIl44,PASSI-<SNSHATf  I44,PASS)-TBUFFI |4A) )/VAR 
IFI.NOT.TYPOUTIGO  TO  62 
MRITEI6,133)PASS,ISNSMAT<I,PASS) ,I-1,B) 

133  FORMAT!*  SENSMAT !  I  , *  ,  12, * )  -  *,BF11.5) 

62  PASS-PASS*! 

GO  TO  132,35, 37, 39, 41, 43, 45, 47, 53), PASS 
C 
C 

C  FINISHED  WITH  AN  ENTIRE  9-PASS  ITERATION 
C 
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53  CCNTINUE 

MRITE(6.134I 

134  FCRMATI//.*  SENSMAT  >  S/l 
CC  66  I  ■UK 

66  hRITE(6t 104) I SNSMAT ( I t J ) t J- 1*B) 

104  FCRMATt 1X.8F11.5) 

CALL  MATINV(SNSMAIt8tl6t 
WR1TE(6.135) 

135  FCRMATI/,'  (SENSMATt««-l  «  •,/) 

DC  67  l>1.8 

67  taRITEI6tl04)(SNSMATII  tJ)tJ-<)«16) 

CC  49  I>1.8 

GRAOl I )>0. 

OC  49  J-ltB 

49  GRAOII )>CRAOII )♦SNSMATIItJ♦8|4ITL0A0IJI-TBUFFIJII 

hRITEI6.136) 

136  FCRMATI/.*  GKAOIENT  VECTOR  ■  './) 
taR|TE(6,104)IGRA0(I ).l-1.8) 

11ER«ITER«1 

MAG>0. 

00  87  l>1.8 

87  MAC>MAG*GRAOI I )**2 
MAG>SQRT(NAG) 

IFIMAG.GT.MAXJBPISTPS1E>MAXJMP/MAG 
IFIMAG.l.C.MAXJMP)SrPS7E«1.0 
OC  50  I>1«11 
CC  50  J>U4 

50  2STRT(  I  )-ZSTRT(  n»ICRAOU)«TCHEB(  I  .  J)  ♦  JJ4CR  A0IJ»4  )  *TCHEBI  I .  J  I  )  * 
1STPS2E 

HRirE(6.137)  STPSZE 

137  FCRMATI/,'  STEPSI2E  MULTIPLIER  •  *.^8.4) 
hRITE(6tl01) 

101  FCRMATI/,'  LOCATION  CAPACITY  2N0  VARIATION',/) 

)<RITEIA,102)IShCHPT(l),CAPII),SECVARII),l«l,NSTEPS) 

102  FCRMATI 1X,F 10.6, 8X,E 15.4, 5X,E 15.4) 

IFIMAG.Lr.lO.E>24MAXJMP)GO  TO  27 
IFI  ITER.LE.ITMAX)  CO  TO  51 

C 

C 

C  CCNE  talTH  THE  ENTIRE  PROBLEM-PUNCH, THEN  READ  ANOTHER. 

C 

C 

52  CCNTINUE 

hRITEl6,145) 

145  FCRMATI//,'  THIS  PROBLEM  IS  FINISHED  ...  WILL  READ  THE  NEXT 
IPRCBLEM  ...  ',//) 

IFI .NOT. PUNCH)  GO  TO  70 
SELFST>. FALSE. 

CCNTIN>. FALSE. 
ltRITEI7,CNTRL) 

HR|TEI7,PARAMS) 

IFIWAVGOE)NRI TEI7,hlOTHl) 

)«RITEI7,ZSTART) 

CC  TO  70 
C 
C 

C  ITERATION  HAS  CONVERGED 

C 

27  NR|TEI6,14e) 

146  FCRMATI//,'  THIS  ITERATION  HAS  CONVERGED...  HILL  READ  THE  NEXT 
IPROHLEM  ...  ',//) 

CC  TO  70 
C 
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c 

c  self-startup  procedure 

72  CONTINUE 
C*CNIN 

CC  13  I-l.Il 
XCLO(I«<i)-l. 

XCLO(It3)-AINAC(2LOAO( I  1) /REAL ( ZLOAOI 1  I  I 
XCL0(It2)-l«/(SQRT2*REAL(ZL0A0(n)) 

XOLOII t l)«XOLO(I t2IP(REAL(ZLOAO(n )**2«A1MA0I ZLOADI I ) 1**21 

if(mavgoe)  go  to  1A 

Tl-C 

GO  TO  13 

lA  Tl-C-l./(L*(2.*WIOTH*FREQ( 111**2) 

13  RCHARI  n-SORT(L/Tl) 

0IST2—LENGTH 

CALL  STEP(OIST2tCMULT) 

OC  IS  I-l.ll 

15  ZSTRTd )-( 1.0*JJ*XNEN( I,3) )/(S0KT2*XNEW( I ,2)) 

GC  TO  73 
C 
C 

C  ERROR  DIAGNOSTICS  FOR  READ 
C 

90  NRITE(6>129I 

129  rCKMATC  ERROR  WHILE  READING  OATASET  /POLY/') 

GC  TO  71 

91  WRITE(6.129) 

125  FCRMATI*  ERROR  WHILE  READING  OATASET  /CNTRL/') 
REAC(6«PARAMStERR>92) 

NERRS>NERRS*1 
IF(NERRS.GE.4)G0  TO  71 
GC  TO  70 

92  NRITE(6.126) 

126  FORMATI'  ERROR  WHILE  READING  OATASET  /PARAMS/') 
NERRS«NERRS*1 

1F(NERRS.6E.A)G0  to  71 
GC  TO  70 

93  WRITE(6.127) 

127  FCRNATC  error  WHILE  READING  OATASET  /ZSTART/') 
NERRS*NERI\S*1 

IF(NERRS.GE.A)G0  to  71 
GC  TO  70 
C 
C 

C  FINAL  EXIT 
C 

71  CONTINUE 
86  WRITE(6.138) 

138  FORMAT!///, •  ••*•••  THATS  ALL  FOLKS  . ,///) 

C 

C 

C 

STOP 

ENO 
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SL9R0UT1NE  0 1 STl ( CKULT ,LOC 1 01  ST .LENGTH ) 

CCHMON  NVAL S.THRESH.L.C, PI ,SCRT2,F0RFLC, SLOPE, CURV,^■ftXSTP. 

IFREQI  ID.XNEWI  ll.Al.PNEWIll.A) 

REAL  LiLENGTH.LCC.HAXSTP 
LCGICAL  FuRFLG 
C  FLAG  TO  SCAN  FCKhAKU 

IFIFORFLGI  FCRSGN-l. 

C  ELSE  SCAN  BACKWARD 

IFI.NOT.  FORFLG)  FCRSGN*-!. 

SLCPE*C. 

CLiRV>0. 

Y5«0. 

Y6«C. 

CC  1  l>lfU 

SLOPE-SLOPE* (XNEWI 1,1 »♦PNEW^I ,1 »-XNfcW( I ,2 ) *nNf  V.  I  1 , 2 ) ) ^FRHQt I ) **2 

1  CLlRV-CLKV*(L*(XNEHI I ,2)*PNEW( I ,  1)*XNEWI I • T)*PNbWI I , 1) >*C* 

IIXNEWI  1 ,  1  )*PNew(  1  ,  3)*XNEWI  I  ,  3 )  *PNEW  I  I  ,  2  )  >  )  *FHC(.' 1 1  )  **3 

SLOPE-fl. •PI ••2*1 ‘SLOPE 
CLRV  — 16*Pl**3*L*SQRT2*CtRV 
C 
C 
C 

IF (NVALS.lt. 3  .OK.  THRESH. EO. 0.0)  GC  TC  12 
IF(CMULT-THRESH)  2,2,3 
C  CPULT  ABOVE  THRESHOLD  REGION 

3  YJ-CHUL T-THRESH 

CO  TO  S 

2  IF(CMULT«THRESH)18,7,7 

C  CPULT  BELOW  THRESH.  REGION 

18  Y7«CMULT*THRESH 

GC  TO  5 

C  CMULT  INSIDE  D'RESH  REGION 

C 
C 
C 

7  Y7«0. 

T3«2.*THRESH/(NVALS-2.) 

Ti«CMULT*INVALS-2.)/(2.*THRESH)*NVALS/2.0-INT(CMULT*(NVALS-2. )/ 
1(2.*THRESH)*NVALS/2.0I 
IFISL0PE*F0RSGM22.23«23 

22  Y5-T2PT3 
GC  TO  2A 

23  Y5«l l.0-T2l*T3 

2A  IF(CURV)2S,26,26 

25  Yt«T2*T3 
GC  TO  27 

26  Ye«( 1.0-T2I*T3 

27  CCNTINUE 
T1-Y6 

C  LINEAR  EXTRAPOLATION 

Y2-F0RSGN*Y5/ABS(SL0PEI 
GC  TO  13 
C 
C 
C 

C  ONLY  TWO  capacity  VALUES 

12  Y7-CMULT 

5  T1-ABS(Y7I 

IF(Y7*SL0PEPF0RSGN)8.10,10 
10  Y2-F0RSGNPMAXSTP 

GC  TO  9 

C  LINEAR  EXTRAPOLATION 

a  Y2  —  Y7/SL0PE 
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9  lF(Y7*CURVI13.11tU 

11  Y1-F0RSGN*MAXSTP 

GC  TO  14 

C  ClAORATIC  EXTRAPOLATION 

13  Yl*-SLOPE/CURV*FORSGN/ABS(CORVI*SORT(SLnPE**2*2.*ARS(CURVI*Tll 

C  PICK  SMALLER  OF  TMO  STEPS 

14  IF(  ABS(  Y1I-ABS(Y2M  15.  15. 16 

16  CIST-Y2 
GC  TO  17 

15  DIST«Y1 

C  TRUNCATE  IF  STEP  IS  TOO  LARGE 

17  IF( ABSIOISTI.GT.MAXSTPI  DIST*F0RSGN*MAXSTP 

C  TRUNCATE  IF  ROUNDOFF  ERRORS  GIVE  WRONG  SIGN 

IF(FORSGN*OIST.LT.  0.01  0IST*0.0 
C  CHECK  FOR  RUNNING  OFF  END  OF  LINE 

IFKLOCtOIST  .GT.  0.0 1  .  AND.  FORFLGI  DIST*-LDC 
C  CHECK  FOR  RUNNING  UFF  BEGINNING  OF  LINE 

IFKLOC^OIST  .LT.  -LENGTH!. and.. not. FORFLGI  DI  ST*-LENGTH-LDC 
C 
C 
C 

RETURN 

EKO 


SUBROUTINE  STEPI0IST2.CMULTI 

CCMMON  NVALS.THRESH.L.C.PI.SURT2.F0RFLG.SL0PE.CURV.MAXSTP. 

IFREOt 11) .XNEWI 11.41 .PNEWI 11,4! .XOLOl 11.41 .POLO! 11.4I.RCHARI 11! 
CCUBLE  PRECISION  TOl, T02, T03. T04. THETA 
LOGICAL  FORFLG 
REAL  LiMAXSTP 
CMULT«0. 

SL0PE>0. 

00  1  I-l.ll 

ThETA>4.«PI«DIST2*L*FREti(I I/RCHARII ! 

TCl-OSINITHETA! 

TC2«0C0S( THETA! 

TC 3*1. 4702 
TC4-1.-TD2 
TS2-RCHAR( I !/S0RT2 
TS1«TS2**2 

TS4*1./(RCHAR( I !*S0RT2! 

TS3»TS4**2 

XKEWU,1!>.5*T03*X0L0(I ,1!4TS1*T04*X0L0(I,2!-TS2*T01*X0L0( 1,3! 
XNEW(I.2!-TS3*T04*X0L0n  ,l!4.5*T03*X0L0n,2!4TS4*T01*X0L0(  1,3! 
XKEW(I,3!-TS4*T01*X0L0(I,1!-TS2*T01*X0L0(I.2!4T02*X0L0(I,3) 

PKEWI I,1!-.5*T03*P0L0( I,1!4TS3*T04*P0L0(I ,2!-TS4*T01*P0L0( 1,3! 
PKEMn,2!-TSl*TO4*POL0(I  ,l!4.5*T03*P0L0n,2!4TS2*T01*P0L0(  1,3! 
PNEWI I ,3!>TS2*T01*P0L0( I , 1 !-TS4*T01*POLO( I , 2 ! 4T02*P0L0( I , 3 ! 
SL0PE>SL0PE4(XNEW( 1,1! •PNEWI 1,1 !-XNEW( 1,2 !*PNEW( 1,2! !*FREO( 1 !**2 
1  CMULT>CMULT42.*S0RT2*P1*FREQI I !*IXNEWII,3!*PNEWI1,2)4 
1XNEWI1,1!*PNEWI1,3!! 

SL0PE»8.*PI**2*L*SL0PE 

C 

C 

c 

RETURN 

ENO 
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SUBROUTINE  SWl  TCH  (CHUL  T  tLOC  tOI  ST  ,  NEk<V  AL  I 

CCMhON  NVALSt  THRE SH , L t C t Rl t SORT 2 « F ORF LG . SLOPE . CUR V , NAX STP, 
IFKECI  111 .XNEhl ll.AI .PNEWI ll.Al.XULOl 1 1 .A  I , PCLDI 1 1 t<>  1 1 RCmAR I  11), 
2NSTEPStOIST2tCF'IKfCHAX.V,AVCOE  .H I DTH ,C  AP I  200 1 ,  SWCHPT  I  200  I 
NAMELIST/FC0VAK/SL0PE,CURV,D1ST,NEWVAL .Tl.rjSTCPi.LOC.C.CMULT 
REAL  LtLOCtMAXSTP 
LCCICAL  WAVCOE.FCRFLC 
IFI THRESH* IN VALS-2 111. 1.2 

1  IFISLOPEI3.3.A 

3  NERVAL>0 
CC  TO  5 

4  NEMVAL>NVALS-1 
CC  TO  •» 

C  MAP  THRESHOLD  REGION  INTO  10.11 

2  Tl-ICMULT/THRESH*!. 1/2. 

NEWVAL>INT(  1.*(NVALS-2I*TU.S*SICNI  1.0.  SLOPE  I  I 
IF INEWVAL.LT.0INEMVAL>0 
IF (NEWVAL.GT.NVALS-1INEWVAL>NVALS-1 

5  NSTEPS-NSTEPS*! 

S*CHPT(NSTEPSI«LCC 

FORFLC-.TRUE. 

0IST2-0. 

C«CM|N«I (CMAX-CMINI/INVALS-ll l•NEWVAL 

CAPINSTEPSI-C 

CURV-0. 

OC  8  I ■  1.11 
IFIWAVCOEI  CO  TO  6 
T2«C 
CC  TO  7 

6  T2«C-l./a*<  2.*MI0TH*FREQI  1 1 1**/! 

7  RCHARI I I>SQRT(L/T2I 

8  CLRV«CURV*(L*(XNEHI|.2I*PNEI)II.3I*XNEWI  l.3l*PNEHn.ll  I^C* 
l(XNEM(I.ll*PNEh(l.3l*XNEWII.3l*PNEWn.2lll*FRECI I  1**3 

CURV«-16.*PI**3*L*S(3RT2*CURV 

C 

C 

C 

IF(THRESH*(NVALS*2I.EQ.0.0I  CO  TO  13 
IP(NEHVAL*(NEWVAL*1->NVALSI  19.13.13 
C  OUTSIDE  THRESHOLD  REGION 

13  IF(SLOPE*CURVI 10.11.11 

10  OIST**2.*SLOPE/CURV 
CC  TO  12 

11  OIST-MAXSTP 
CC  TO  12 

C  INSIDE  THRESHOLD  REGION 

9  0IST>2.0*THRESH/i INVALS-2I*ABSI SLOPE  I  I 

12  IFIABSIOISTI.CT.MAXSTP)  OIST>MAXSTP 
IFI IL0C*01STI.CT.  0.0)  OIST>-LOC 

DO  14  I>1.11 
00  14  J*1.4 
XCL0I1.JI>XNEW(I.JI 

14  PCLOl I.J)>PNEHII.J) 
IPICAPINSTEPSI.EQ.CAPINSTEPS-lllWRITEIN.FUUVARI 
RETURN 

END 


SLBRQUTINE  VALUE  I CMULT .NEWVAL ) 
CCMHON  NVAL St  thresh 
lFITHRESHAINVALS-21 )6.6.7 
Tl-ICHULT/THRESH«l.l/2. 

GC  TO  8 

Tl-SlCNI2.0fCMULT) 

lF(Tni,l,2 

NEMVAL-0 

GC  TO  3 

IFIT1-I.I4.4.S 

NENVAL>NVALS-l 

GC  TO  3 

NEhVAL-INTI  (hVALS-2>ATUl.l 
CONTINUE 


RETURN 

EKO 


103 


r»  f  ^  r>  r> 


SUBROUTINE  MAT INV ( M .ROMS tCOLSI 
INTEGER  ROMS. COLS 
INTEGER  Cllll 
REAL  MIROMS.COLSI.T2(22l 

C  ENTER  MITH  MII.JI.I*!  TO  ROMS*  J«1  TO  ROMS 

00  1  1-1. ROMS 

C  ADJOIN  lOENTITV  MATRIX 

Cll l-l 
IRPl-ROMS*! 

OC  1  J-IRPl.COLS 
IFIJ-l-ROMS)2.3.2 

2  MII.JI-0. 

GO  TO  1 

3  MII.J)-1. 

1  CONTINUE 


GAUSS  REDUCTION  ON  EACH  ROM 
00  10  K-l.ROMS 
C  PINO  MAX  ABS  ELEMENT  AVAILABLE  FOR  REDUCTION 

ABSMAX-O. 

00  S  1-1. ROMS 
00  5  J-1*R0MS 
IFICIinS.S*6 
6  IFICIJIIS.S.r 

T  IFIABSIMII. JII-ABSMAXIS.S.e 

8  ABSMAX-ABSIMII.JII 
IMAX-I 

JMAX-J 

5  CONTINUE 

C 
C 
C 

CIJMAXI-0 
Tl-MI IMAX.JMAXI 

C  SNAP  RONS-MAX  ABS.  ELEMENT  ONTO  DIAGONAL 

00  9  J-ltCOLS 
T2IJI-MI IMAX.JI 
Ml IMAX.JI-MIJMAX.JI 

9  Ml JMAX.J)-T2IJ)/T1 

C  GAUSS  REDUCTION  ON  COLUMN  JMAX 

00  10  1-1. RONS 
lFll-JMArill.10.11 
11  Tl-MI I.JUAXI 

OC  10  J-l.COLS 
MII.JI-MII. JI-T1*M( JMAX.JI 

10  CCNTINUE 
C 

C 

C 

RETURN 

C  LEAVE  MITH  MII.JI.I-1  TO  ROMS.  J-RONS^l  TO  COLS 
C 

c 

c 

ENO 
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within  0.2  percent  of  that  obtained  using  a  continuum  of  dielectric  constants  I^i^lN’^MAJU 
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